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FOREWORD 


The  investigation  reported  herein  is  part  of  the  research  project  at 
The  Ohio  State  University,  Columbus,  Ohio  supported  by  the  Air  Force  Of¬ 
fice  of  Scientific  Research  Grant  83-00-55.  Lt.  Col.  John  J.  Allen  was 
the  Program  Manager  at  the  commencement  of  the  Project.  Lt.  Col.  Law¬ 
rence  D.  Hokanson  was  the  Program  Manager  from  July  1,1983.  The  Research 
project  was  started  in  February  1,  1983  and  is  continuing.  The  present 
report  documents  part  of  the  work  done  up  to  January  31,  1984.  At  The 
Ohio  State  University,  the  Project  is  being  supervised  by  Dr.  Ranbir  S. 
Sandhu,  Professor,  Department  of  Civil  Engineering.  The  computer  program 
modification  reported  herein  were  carried  out  by  Dr.  Baher  L.  Aboustit, 
Post-doctoral  Research  Associate.  The  documentation  of  this  report  was 
done  by  Graduate  Research  Associates  Soon-jo  Hong  and  Mahantesh  S.  Hire- 
math.  The  Instruction  and  Research  Computer  Center  at  The  Ohio  State 
University  provided  the  computational  facilities. 


ABSTRACT 


A  computer  program  was  developed  for  evaluation  of  finite  ele¬ 
ment  models  for  soil  consolidation  and  study  of  dynamic  response  of  flu¬ 
id-saturated  soils.  One-  and  two-dimensional  consolidation  problems  were 
analysed  using  different  finite  elements.  Transient  response  of  saturat¬ 
ed  porous  elastic  media  for  dynamic  as  well  as  quasi-static  problems  was 
studied  Results  were  compared  with  the  numerical  and  analytical  solu¬ 
tions  available. 
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SECTION  I 


INTRODUCTION 

Numerical  performance  of  Ghaboussi  and  Wilson's  [1]  isoparame¬ 
tric  bilinear  quadri 1 ateral  element  was  compared  with  that  of  Sandhu 
[2,3,4]  for  analysis  of  quasi-static  flow  of  an  incompressible  fluid 
through  a  linear  elastic  saturated  soil  in  reference  [5].  Both  these  el¬ 
ements  were  used  to  solve  Terzaghi's  problem  of  one-dimensional  consoli¬ 
dation  and  Gibson's  [6]  analysis  of  two-dimenSional  plain  strain  consol¬ 
idation.  The  results  from  the  computer  program  developed  for  this 
purpose  were  compared  with  the  well-known  theoretical  solutions. 

Reference  [7]  used  one-dimensional  linear  element  and  4-node  and 
8-node  isoparametric  two-dimensional  elements  for  numerical  solution  of 
several  one-dimensional  problems  for  ouasi-static  as  well  as  dynamic 
loading  conditions.  Here  again,  the  results  of  the  computer  analysis 
were  checked  with  the  exact  solutions  available  in  literature. 

This  report  forms  the  continuation  of  the  investigation  carried 
out  in  [5]  and  [7].  The  computer  program  developed  for  implementaion  of 
the  theory  presented  in  [5]  and  [7]  is  described  here.  The  program  is 
capable  of  handling  one-  and  two-dimensional  problems  for  steady-state 
and  transient  solutions  for  soil  consolidation.  Quasi-static  and  dynamic 
loading  conditions  have  been  included.  Compressible  as  well  as  incom- 


pressible  fluids  are  permitted.  Time  domain  integration  is  carried  out 
by  using  Wilson's  [1]  (3-^-0  method.  However,  the  program  is  limited  for 
linear  analysis  of  soil  consolidation  i.e.  small  deformations,  linear 
elastic  soil  behavior,  constant  permeability  and  satuzated  soils. 

The  computer  program  described,  is  writen  in  FORTRAN  IV  language, 
with  which  engineers  are  well  versed.  The  program  is  modular  and  versa¬ 
tile  enough  to  classify  as  a  research  tool,  in  which  'modules'  can  be 
changed  or  added  as  desired.  In  order  to  have  a  computer  program,  which 
can  solve  many  types  of  finite  element  problems,  a  macro  programming 
language  is  introduced.  The  macro  programming  language  ia  associated 
with  a  set  of  compact  subprograms ,  each  desired  to  compute  one  or  at 
most  a  few  basic  steps  in  a  finite  element  solution  process.  This  con¬ 
cept  also  permits  inclusion  of  different  solution  algorithms.  Dynamic 
storage  allocation  has  been  used  to  economise  on  storage.  Thus,  each  ar¬ 
ray  is  variably  dimensioned  by  using  a  set  of  pointers  established  in 
the  control  program. 

Section  II  discusses  the  program  capability  including  types  of 
problems  analysed,  elements  used,  field  equations  adopted  and  solution 
process  i ncorporated.  Section  III  gives  program  structure  with  flow 
charts  for  individual  problems,  nesting  of  subroutines  along  with  brief 
discussion  on  each  subroutine.  Section  IV  includes  the  instructions  for 
program  usage,  particularly  data  input  module.  A  complete  listing  of  the 
program  is  placed  in  the  Appendix. 


SECTION  II 


PROGRAM  NARRATIVE 

2.1  PURPOSE  OF  THE  COMPUTER  PROGRAM 

The  computer  program  obtains  numerical  solution  using  finite  ele¬ 
ment  method  for  one-  and  two-dimensional  consolidation  problems.  The 
details  are  available  in  [5],  Displacements  and  pore  water  pressure  are 
the  field  variables  used  in  these  problems.  The  program  also  analyses 
the  dynamic  response  of  fluid  saturated  porous  media.  Quasi-static  anal¬ 
ysis  is  also  included.  Solid  displacements  and  relative  displacement  of 
fluid  with  respect  to  soil  skeleton  are  used  as  unknown  quantities  in 
these  problems.  The  time  domain  integration  is  carried  out  by 
method  [7], 

2.2  METHOD  OF  ANALYSIS 


The  finite  element  is  used  for  obtaining  numerical  solutions.  The 
method  is  wel 1 -documented  in  literature  and  details  of  these  are  not  in¬ 
cluded  in  this  report.  Variational  principle  leading  to  matrix  fomula- 
tion  used  in  case  of  consolidation  analysis  is  represented  in  [5],  while 
Galerkin  finite  element  method  used  in  response  study  of  saturated  soils 
is  presented  in  detail  in  [7],  The  theoretical  background  for  time  do¬ 
main  integration  using  Wilson's  [1]  8-^-0  method  is  also  given  in  [7]. 
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2.3  FINITE  ELEMENTS  USED 


For  soil  consolidation  analysis  [5],  two  types  of  elements  are 
used  viz.  Ghaboussi  and  Wilson's  [1]  6-4  element  and  Sandhu's  [3]  8-4 

composite  element.  Ghaboussi' s  element  is  more  like  4-4  isoparametric 
element,  but  uses  additional  'local'  nodes  that  give  it  a  character  of  a 
'higher  order  scheme*.  Sandhu's  element  uses  8-point  biquadratic  inter¬ 
polation  for  dispalcements  and  4-point  isoparametric  quadrilateral  for 
fluid  pressure.  In  response  analysis  of  fluid-saturated  soil  in  [7], 
three  types  of  elements  have  been  used.  These  are  one-dimensional  linear 
element,  4-node  i soparametric  element  and  8-node  isoparametric  element. 
Details  about  finite  element  formulations  for  these  two  analyses  and  el¬ 
ement  development  is  given  in  [5]  and  [7]. 

2.4  ASSUMTIONS  MADE  BY  THE  PROGRAMMER 

In  consolidation  analysis  [5],  the  fluid  was  assumed  to  be  incom¬ 
pressible  in  solving  both  the  one-  and  two-dimensional  problems.  For 
other  material  properties,  refer  to  [5].  In  2-dimensional  analysis, 
Poisson's  ratio  was  set  equal  to  zero.  Also,  the  finite  domain  was  mod¬ 
eled  by  a  finite  one  in  the  numerical  schemes.  Following  alternative 
boundary  conditions  were  specified  at  the  boundary,  where  half  space  was 
cut  off. 

1.  Pore  fluid  pressure  and  horizontal  displacements  prescribed. 

2.  Fluid  flux  and  horizontal  displacements  precribed. 
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3.  Pore  presure  and  traction  prescribed. 

In  temporal  discretization,  piecewise  linear  variation  of  displacement 
and  pore  fluid  pressure  was  assumed. 

In  the  report  on  response  analysis  of  fl  ui d-saturated  porous  medi¬ 
um  [7],  following  problems  were  solved  using  the  computer  program  : 

1.  Quasi-static  soil  consolidation. 

2.  Dynamic  response  of  an  elastic  layer  of  single  material. 

3.  Response  of  fluid-saturated  soil  layer. 

All  the  problems  deal  with  the  compressional  wave  propagation  in  an  ini¬ 
tially  undisturbed,  homogeneous,  isotropic,  elastic,  porous  or  non-po- 
rous  system  due  to  dynamic  loading.  Solutions  for  non-porous  media  were 
accomplished  by  prescribing  relative  displacement  between  soil  and  fluid 
to  be  zero.  Quasi-static  consolidation  problems  were  solved  by  setting 
density  and  damping  to  zero.  For  Ghaboussi's  [1]  problem,  damping  was 
assumed  to  be  linear  combination  of  stiffness  and  mass  matrix  as  given 
in  [1],  The  time  domain  integration  scheme  using  Wilson's  [1] 
method  utilized  various  values  of  0  and  0  for  stability  analysis.  In 
analysis  of  Garg's  [7]  problems  with  'strong'  and  'weak'  coupling  dif¬ 
ferent  values  of  permeability  and  time-steps  were  used.  Further  theoret¬ 
ical  details  are  given  in  [7]. 


2.5  FIELD  EQUATIONS  OF  ANALYSES 


2.5.1  Equations  Governing  Linear  Elastic  Soil  Consolidation 

Assuming  pore  water  to  be  compressible,  the  equations  of  force  equi¬ 
librium  of  elementary  volumes  and  mass  continuity,  over  the  spatial  re¬ 
gion  of  interest  R,  may  be  written  in  standard  indicial  notation  as 


^•Eklijuk,l^,i  +  +  f  j  0 

(1) 

*-Kij^,i  +^2fi^J  +  uj,j  = 

(2) 

where  u^  f . t  E^^-j,  K^j,  denote  the  cartesian  components,  respectively, 
the  displacement  vector,  the  body  force  vector  per  unit  mass,  the  isoth¬ 
ermal  elasticity  tensor  and  the  permeability  tensor.  ^  is  the  mass  den¬ 
sity  of  the  saturated  soil  and  ^at  of  water.  If  is  the  pore  water 
pressure,  oCis  the  solid  compressibil i ty  and  M  is  a  measure  of  fluid 
compressibility.  With  these  field  equations  we  associate  the  following 
boundary  conditions; 

a 

u,-  on  S^i 

^ijnj  =  ti  on  s2i 
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Here,  t- t  are  components  of  the  traction  and  fluid  flux  vectors  asso¬ 
ciated  with  surfaces  embedded  in  the  closure  of  R.  T—  are  components  of 
the  total  stress  tensor.  $2,  are  complementary  subsets  of  the  bound¬ 

ary  of  the  spatial  region  of  interest  and  so  are  S^,  S4.  Even  though  the 
equations  given  above  apply  to  compressible  fluids,  the  applications  re¬ 
ported  herein  assumed  incompressible  fluid  i.e.  M  ->  oO  .  The  initial 
conditions  for  the  problem  are 

ui (xj ,  0)  =  Ui (0)  on  R 
If  ( x  . ,  0)  =  (0)  on  R 

2.5.2  Equations  Governing  Dynamics  of  Fluid-Saturated  Media 

Biot's  [8]  equations  of  motion  for  an  elastic  porous  medium  saturated 
with  a  compressible  fluid  may  be  written  in  standard  indicial  notation 
as; 


^i  jkluk,l  +  *  M^uk,k  +  wk,k^,j  +^fi  =  ^  u’i  +  "  -  ^  2  W,*  (4) 

i  «  <1  1  • 


[M(0(uk,k + 


“k,k»,f  *  -j-  Vi  ■  V.- +  72 ?2"'  *  r  <S) 


where  u^  f  W]- ,  f^ ,  E^j^i  denote  the  cartesian  components,  respectively, 
of  the  solid  displacement  vector,  the  relative  fluid  displacement  vector, 

the  body  force  vector  per  unit  mass  and  the  isothermal  elasticity 
tensor.  8  is  mass  density  of  the  saturated  soil  and  8-,  that  of  water  per 
unit  bulk  volume,  f,  K,  o(,  M  are,  respectively,  the  porosity,  the  perme- 


ability,  the  solid  compressibility  and  the  fluid  compressibility.  The 
superposed  dot  implies  a  time  derivative.  All  the  functions  are  defined 
over  the  cartesian  product  R  x  [0,#>)  where  R  is  the  spatial  region  of 
interest  and  [0,00)  is  the  positive  interval  of  time.  With  these  field 
equations,  we  associate  the  following  boundary  conditions. 


on  Sj-j 

(6) 

ijkluk,l  + 

o(f[6ij)nj  =  tj  on  S2l- 

(7) 

A 

7C 

+ 

at 

II 

TT  ( t )  on  S3 

(8) 

on  S4 

(9) 

where  S^. ,  s2l-  are  complementary  subsets  of  the  boundary  S  of  the 
spatial  region  of  interest  and  so  are  S3,  s4.  The  initial  conditions  for 
the  problem  are  given  by: 
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2.6  SOLUTION  PROCESS 


2.6.1  Consol idation  Analysis 

Once  the  finite  element  discretization  is  done,  the  solution  pro¬ 
cess  is  initiated  by  setting  up  the  'spatial  stiffness  matrix',  'spatial 
flow  matrix'  and  'coupling  matrix'  for  each  element.  For  further  details 
refer  to  [5],  The  element  load  vectors  are  also  set  up  as  given  in  [5], 
After  assembly  is  done  as  usual,  geometric  boundary  conditions  are  ap¬ 
plied  and  solutions  for  displacements  and  pore  fluid  pressures  are  ob¬ 
tained  for  the  load  under  consideration.  The  load  vector  is  updated  for 
subsequent  time  intervals  and  the  process  is  continued  till  all  time-in¬ 
crements  are  completed. 

2.6.2  Dynamic  Response  Analysis 

Following  the  finite  element  formulation  presented  in  [7],  the 
first  step  in  dynamic  response  analysis  is  setting  up  of  element  stiff¬ 
ness,  mass  and  damping  matrices  and  load  vectors.  The  standard  assembly 
process  follows.  The  time  domain  integration  is  used  to  solve  for  solid 
displacements  and  relative  fluid  displacements,  which  are  the  unknown 
quantities  for  this  analysis.  This  is  used  for  calculation  of  stresses. 
The  process  is  repeated  for  successive  time  intervals,  while  values  of 
load  vector  are  updated  at  the  end  of  each  time  interval. 


SECTION  III 


PROGRAM  STRUCTURE 

3.1  INTRODUCTION 

The  computer  program  developed  for  soil  consolidation  and  dynamic 
response  analysis  is  general  in  nature  and  modular  in  structure,  each 
module  is  designated  to  carry  out  one  or  at  most  a  few  basic  steps  in 
finite  element  solution  process.  The  use  of  macro  programming  language 
is  ideally  suited  for  a  such  a  comprehensive  program,  which  is  capable 
of  solving  many  types  of  finite  element  problems.  Individual  modules  can 
be  easily  modified  or  added  to  further  the  capabilities  of  this  program. 
The  use  of  dynamic  storage  allocation  permits  economic  use  of  storage 
capacity.  The  program  is  user  friendly  and  easy  to  follow  once  macro 
programming  concept  is  understood. 

3.2  FLOW  CHARTS 

The  flow  charts  for  soil  consolidation  and  dynamic  response  analy¬ 
sis  are  shown  in  figures  (1)  and  (2). 
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r  > 

START 


Read  and  print  data 


Compute  the  elementsti f fness ,  mass,  damping 
matrices  and  load  vectors  and  save  it  assemble 
to  get  the  global  [K]  and  [R] 


Obtain  global  mass  matrix  and  initial 
acceleration 


Multiply  element  damping  and  mass  by  the 

appropriate  coefficients  of  P-v-e  method 

* 

and  add  it  [_KJ  to  obtain  [K  ]  and  apply 
geometric  boundary  condition 


Triangulize  [K] 


Update  (RiJ 


Solve  for  {00^  {v**}1  and  compute  stresses 


. _ l _ „ 

<  Are  all  time  increment  completed  > 


Is  t  changed 


Figure  2:  Flow  Chart  For  Dynamic  Response  Analysis 


k 

p 


3.3  NESTING  OF  SUBROUTINE 


The  various  subroutines  used  and  their  structural  relationship  is 
shown  in  figure  (3). 

3.4  SUBROUTINES 


Following  paragraphs  outline  the  operations  performed  in  each  pro¬ 
gram  uni  t. 

3.4.1  MAIN 

In  this  unit,  the  program  capacity  is  set  up.  It  calls  the  subrou¬ 
tine  ACNTRL,  which  is  the  principal  control  unit  of  the  entire  program. 

3.4.2  Subroutine  ACNTRL 

This  subroutine  is  called  by  the  unit  MAIN.  It  forms  the  central 
core  unit  of  the  program  and  calls  several  other  subroutines  to  carry 
out  specific  operations.  It  also  reads  in  the  control  data  and  sets  the 
pointers  for  dynamic  storage  allocation. 

3.4.3  Subroutine  ACTCOL 

This  subroutine  is  called  by  ACNTRL  and  is  in  active  column  pro¬ 
file  symmetric  equation  solver.  The  columns  above  the  diagonal  or  the 
rows  below  the  diagonal  are  stored  in  a  single  subscript  array  and  a 
pointer  is  used  to  locate  the  diagonal  elements. 


a 

SJ 
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3.4.4  Subroutine  ADDSTF 

This  subroutine,  called  by  subroutine  FORMEL,  assembles  the  global 
arrays  by  direct  stiffness  method.  This  assembly  is  for  stiffness,  mass 
and  forces.  Logical  variables  are  set  to  define  which  array  or  arrays 
are  to  be  assembled. 

3.4.5  Subroutine  BANADD 

Subroutine  FORMEL  calls  this  unit  to  assemble  the  element  arrays 
(EK)  and  (EF)  to  form  the  global  arrays  viz.  ( GK )  and  (GF).  The  global 
arrays  are  stored  in  compacted  form. 

3.4.6  Subroutine  BANSOL 

This  subroutine  is  called  by  ACNTRL  and  CONSOL  subroutines.  It 
triangul arises  a  banded  and  symmetzic  coefficient  matrix  and  then  solves 
the  system  of  linear  equations. 

3.4.7  Subroutine  BCCODE 

This  subroutine  is  called  by  ACNTRL  and  reads  and  generates  bound¬ 
ary  constraints  code  for  each  node. 

3.4.8  Subroutine  BNCOND 

This  subroutine  is  called  by  FORMEL,  to  modify  force  vector  ac¬ 


cording  to  element  boundary  conditions. 


3*4.9  Subroutine  CPTABL 

This  subroutine  is  called  by  ACNTRL  to  estimate  time  required  for 
setting  up  stiffness  matrix  for  each  element. 

3.4.10  Subroutine  CONSOL 

This  is  prime  control  unit  for  the  consolidation  analysis.  It  is 
called  by  ACNTRL  and  sets  up  the  finite  element  scheme  by  calling  sever¬ 
al  other  subroutines. 

3.4.11  Subroutine  DOT 

This  function  subprogram  is  used  to  carry  out  the  vector  dot  prod¬ 
uct  and  is  used  by  subroutines  ACTCOL. 

3.4.12  Subroutine  DYNAMC 

This  is  a  prime  unit  for  dynamic  analysis  part  of  the  program  and 
calls  a  series  of  subroutines  to  complete  the  scheme  of  dynamic  analy¬ 
sis. 

3.4.13  Subroutine  ELDATA 

This  subroutine  is  called  by  ACNTRL  to  read  and/or  generate  ele- 


3.4.14  Subroutine  ELMLIB 


This  subroutine  is  called  by  FORMEL  and  is  essentially  a  element 
library  generation  unit.  It  calls  in  element  subroutines  marked  for  par¬ 
ticular  problems. 


3.4.15  Subroutine  ELMT01 


This  subroutine  is  for  Ghaboussi's  element  called  as  6-4  quadri¬ 
lateral  element  used  in  one-dimensional  consolidation  problem. 


3.4.16  Subroutine  ELMT02 


This  element  subroutine  sets  up  the  stiffness,  mass  and  damping 
matrices  in  the  dynamic  analysis  problem.  It  also  evaluates  the  stresses 
at  the  center  of  each  element.  It  is  used  for  one-dimensional  dynamic 
response  analysis  problems. 


3.4.17  Subrouti ne  ELMT03 


This  subroutine  is  similar  to  ELMT10  but  uses  8-noded  quadrilater¬ 
al  element.  It  is  used  for  plane  strain  dynamic  response  analysis. 


3.4.18  Subroutine  FORMEL 


This  subroutine  called  by  ACMTRL  and  loops  over  each  element  to 
set  up  element  matrices  and  later  assemble  them  into  global  arrays. 
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3.4.19  Subroutine  GENVEC 

This  subroutine  called  by  ACNTRL,  formulates  the  nodal  geometry 
data  either  by  reading  it  or  generating  it.  Linear  interpol ation  is  used 
for  this  purpose. 

3.4.20  Subroutine  GEOMBC 

This  subroutine  called  by  TMSTP,  STEPS  and  INACL  and  involves  ap¬ 
plying  kinematic  constraints  for  each  degree  of  freedom. 

3.4.21  Subroutine  GABOSI 

This  is  a  part  of  the  scheme  to  set  up  Ghaboussi's  element  and  is 
called  by  ELMT09  for  tthis  purpose. 

3.4.22  Subroutine  INACl 

This  subroutine  is  called  by  DYNAMC  and  assemles  mass  matrix  and 
solves  the  dynamic  equation  of  motion  for  initial  values  of  accelera¬ 
tion. 

3.4.23  Subroutine  IfJCRS 

This  subroutine  is  called  by  the  unit  DYMAMC  to  set  up  Wilson's 
o-o-O  method  and  updates  the  load  vectors  for  successive  cycles. 


3.4.24  Subroutine  MATLIN 


This  subroutine  is  called  by  ACNTRL  to  read  in  all  material  data 
and  set  up  strain-stress  relations  for  isothermal  material. 

3.4.26  Subrouti ne  MATRLD 

This  subroutine  is  called  by  MATLIN  for  setting  up  material  prop¬ 
erties  for  stress  analysis  problem. 

3.4.26  Subroutine  MELAW 

This  subroutine  is  called  by  MATLIN  for  formulating  stress  strain 
rel ationshi ps  for  isothermal  element. 

3.4.27  Subrouti ne  PCOMP 

This  logical  function  subprogram  is  called  by  ACNTRL,  GENVEC, 
PRTDIS  to  check  for  overflow  during  computations. 

3.4.28  Subroutine  PLOAD 

This  subroutine  called  by  ACNTRL  formulates  the  load  vector. 

3.4.29  Subroutine  PROFIL 

This  unit  is  called  by  ACNTRL  to  compute  the  profile  for  the  glob¬ 
al  arrays  so  that  they  could  be  stored  into  compact  form.  The  equation 
numbers  are  also  set  up. 
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3.4.30  Subroutine  PRTDIS 

This  subroutine  is  called  by  ACNTRL  and  CONSOL  to  print  out  the 
nodal  values  of  the  unknown  quantities  for  each  problem. 

3.4.31  Subrouti ne  PVALUE 

This  subroutine,  called  by  ELMTQ1  and  ELMT03  to  calculate  the  val¬ 
ues  of  principal  strains  and  stresses. 

» 

3.4.32  Subroutine  SET MPT 

This  subroutine,  called  by  ACNTRL,  checks  the  storage  capacity. 

3.4.33  Subroutines  SHAPE2 ,  SHAPE22 ,  SHAPEX 

i 

These  subroutines  are  meant  for  setting  up  the  shape  functions  for 
various  elements,  calculating  the  derivatives  and  transformation  to  nat- 
j  ural  coordinates. 

3.4.34  Subroutines  TMINC 

1  This  subroutine  forms  the  part  of  the  dynamic  analysis  problem  and 

updates  the  load  vectors  for  successive  time  increments. 


L 
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3.4.35  Subroutine  TMSTP 


This  subroutine  does  a  part  of  setting  up  of  stiffness  matrices  in 


SECTION  IV 


COMPUTER  PROGRAM  USAGE 

4.1  PROGRAM  CAPABILITY 

The  code  described  in  this  report  is  a  modified  version  of  program 
"LAMP"  originally  developed  by  the  author  in  the  course  of  his  doctoral 
research.  The  code  was  used  in  [5,7]  to  predict  the  quasi-static  re¬ 
sponse  and  dynamic  response  of  saturated  elastic  media  to  any  boundary 
conditions  or  application  of  external  forces.  The  code  can  handle  me¬ 
chanical  as  well  as  hydraulic  anisotrophy,  and  fluid  compressibility. 
The  program  allows  for  a  choice  of  several  element  types.  These  include 
the  8-8,  the  8-4  and  6-4  elements  for  quasi-static  consolidation  (Figure 
5),  and  the  1-1,  the  4-4  and  the  8-8  elements  for  dynamic  response  anal¬ 
ysis  (Figure  6).  Plane  strain,  plane  stress  as  well  as  axisymmetric  cas¬ 
es  can  be  handled.  Linear  elastic  as  well  as  elasto-plastic  behavior  is 
allowed  for.  The  program  will  solve  one-  and  two-dimensional  mass  and 
heat  transfer,  convective  transport  problems  as  well  as  coupled  deforma¬ 
tion  and  mass  transport.  The  example  problems  stated  in  [5,7]  are  a  spe¬ 
cial  case  of  this  latter  class. 
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4.2  PROGRAM  LIMITATIONS 

For  analysis  of  soil  consolidation,  the  program  capability  is  limited 
to  the  case  of  small  deformations,  linear  elastic  soil  bahavior,  con¬ 
stant  permeability  and  continued  saturation.  The  program  assumes  the  ap¬ 
plied  boundary  conditions  and  loads  constant  throughout. 

4.3  INSTRUCTIONS  FOR  USE 

The  program  is  basically  composed  of  two  modules,  namely  the  data  in¬ 
put  module  and  the  solution  and  output  module.  They  are  controlled  by  a 
set  of  easily  recognized  macro  commands.  This  feature  allows  a  user  to 
program  input  data  mode  and  the  corresponding  flow  of  solution  process. 
Followings  are  the  sequential  order  of  input  data  cards  in  which  dis¬ 
tinctions  are  made  between  quasi-static  consolidation  and  dynamic  re¬ 
sponse  analysis  when  needed. 

4.3.1  Title  Card  (20A4) 


Columns 

Vari able 

Description 

1-4 

TITL  (1) 

Must  contain  LAMP  to  start  the  program. 

5-80 

(TITL  (I), 

Alpha  numeric  title  of  the  problem. 

1=2,20) 

Control 

Data  Card-NAMELIST  Input  with  the 

Name  CONTROL 

Variable 

Descri  ption 

NUMNP 

Total 

number  of  nodal  points 

2 


HU  MEL 


NUMMAT 


NSTEPS 


Total  number  of  elements 

Number  of  different  material  sets 

Number  of  nodes  per  element 

Spatial  dimension 

Number  of  unknowns  per  nodes 

Number  of  pressure  nodes 

Number  of  time  steps 


Note  1;  (a)  If  NAMELIST  input  is  to  be  continued  on  to  the  next  cards,  a 
blank  on  the  first  column  serves  as  a  signal  for  continua- 


(b)  For  quasi-static  consolidation  and  dynamic  response  analysis 
the  following  elements  can  be  selected. 


For  Quasi-static  Consolidation  Analysis 

The  problem  is  solved  only  for  plane  strain  state.  The 
noal  values  are  (ux,  Uy.TT),  i.e.  NDF=3  and  three  types 
of  elements  can  be  used  independently  (Figure  5). 

The  8-4  element  in  which  NEN=8,  NPN-4 
The  8-8  element  in  which  NEN=8,  NPN=8 
The  6-4  element  in  which  NEN=4,  NPN=4 


For  Dynamic  Response  Analysi s 

The  problem  is  solved  for  one-dimensional  and  plane  strain. 
For  one-dimensional  analysis  (NDM=1),  the  1-1  element  where 
the  linear  i nterpol ation  is  employed  is  used.  The  nodal 
values  are  (uv 


wv),  i.e.  NDF=2,  and  NEN=2.  For  plane 


strain  analysis  (NDM=2),  the  nodal  values  are  (u  u..,w¥,w  ) 

x  y  x  y 

i.e.  NDF=4,  and  two  types  of  elements  can  be  selected  in¬ 
dependently  (Figure  6). 

The  4-4  element  which  is  defined  dy  NEN=4 
The  8-8  element  which  is  defined  by  NEN=8 

Example;  the  below  is  an  example  of  NAMELIST  input  for  the  Terzaghi's 
consolidation  problem. 

b&C0NTR0LbNUMNP=48 ,NUMEL=9 ,NUMMAT=1 ,NDF=3 ,NEN=8 ,NPN=4,NSTEPS=5 

4.3.3  Node  Geometry  Data 

4.3.3. 1  Macro-Control  Card 

The  macro  instruction  N0DE(A4)  which  causes  node  input  mode  must  be 
supplied  at  the  first  card. 

4. 3. 3. 2  Coordinate  Data  Card  (215 ,6F10 .0 ) 

Repeat  as  many  times  as  necessary  after  NODE  card. 

Columns  Variable  Description 

1-5  N  Node  number 

6-10  NG  Generation  increment  (see  Note  2) 

11  -  (X( I ,N ) ,1=1 ,NDM)  I-th  coordinate  value  of  node  N 

{ FACTOR (I),I=1,NDM)  Spacing  increment  factor  (see  Note  2) 

Terminate  with  a  blank  card. 


Note  2:  Nodal  coordinates  can  be  generated  based  on  input  values  on  two 
successive  cards.  The  node  number  is  computed  sequentially  (i.e. 
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N  ,N+NG ,N+2NG ,  — .  etc.)  until  the  value  of  N  on  the  second  card 
is  reached.  NG=0  or  blank  causes  no  generation.  Nodal  spacing  in 
I-th  coordinate  of  generated  nodes  will  be  increased  (decreased) 
successively  by  the  F ACTOR ( I )  prescribed  on  the  first  card  of 
the  pair  if  positive  (negative).  FACTOR ( I ) =0  will  result  in 
equal  spacing  in  I-direction.  Nodal  number  N  needs  not  to  be  in 
in  order.  MG  can  be  negative. 

4.3.4  Element  Data 
4.3.4. 1  Macro  Control  Card 

The  macro  instruction  EL£M(A4)  at  the  first  card  causes  element  con¬ 
nectivity  input  mode. 


4. 3. 4. 2  Element  Data  Card  (3I5.NENI5) 

Repeat  as  many  times  as  necessary  after  ELEH  card. 


Columns  Variable 

1-5  L 

6-10  LK  —  I X { NEN 1 , L  ) 

11-15  LX 

16-  ( IX ( J  ,L  )  ,J  =  1  ,NEN ) 


Description 
Element  number 

Material  set  number  ( <NUMMAT ) 
Generation  increment  of  node 
J-th  nodal  number 


Note  3:  Elements  must  be  in  order  and  the  nodal  numbers  must  be  read  in 
counter  clockwise.  If  element  cards  are  omitted,  the  element 
data  will  be  generated  from  the  previous  element  with  the  same 
material  number  and  all  nodes  are  incremented  by  LX  on  the 
previous  element.  Generation  up  to  the  NLImel  occurs  when  a  blank 


card  is  encountered.  No  blank  card  is  needed  if  the  last  data  is 
gi ven. 

Note  4:  The  consolidation  elements  in  Note  1  can  be  generated  by  the 
common  shape  function  routines,  SHAPE2 ,  SHAPE22  and  SHAPEX  ac¬ 
cording  to  number  of  nodes/elements. 

4.3.5  Boundary  Restrai nt  Code  Data 

4. 3. 5.1  Macro  Control  Card 

The  macro  instruction  BOUN ( A4 )  at  the  first  card  causes  boundary  re¬ 
straint  data  input  mode. 

4. 3. 5. 2  Boundary  restrain  code  data  cards  (615) 

Repeat  as  many  times  as  necessary  after  BOUN  card. 

Columns  Variable  Description 

I- 5  N  Node  number 

6-10  NX  Generation  increment  of  node 

II- 10+5*NDF  (10(1),  1=1, NDF)  Restraint  code  in  I-direction 

Terminate  with  a  blank  card. 

Note  5:  Nodal  numbers  are  incremented  by  NX  up  to,  but  not  including  ,  N 
on  the  next  card.  NX=0  causes  no  generation.  I D ( I ) = 1  if  I-th 
degree  of  freedom  is  prescribed,  I D ( I ) =0  if  not.  Only  nonzero 
codes  need  to  be  specified. 
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4.3.6  Nodal  Boundary  Value  Data 

4.3.6. 1  Macro  Control  Card 

The  macro  instruction  3VAL(A4)  at  the  first  card  causes  boundary  data 
i nput  mode. 


4. 3. 6. 2  Macro  Nonzero  Boundary  Data  Card  (215, 7F10 .0 ) 

Repeat  as  many  times  as  there  are  boundary  values. 


Columns  Variable 

I- 5  N 

6-10  NG 

II- 20+10 *NDF  (F( I ,N ) , I =1 ,NDF ) 

( F  ACT  OR ( I ) , I  =  1 , NDF ) 


Descri ption 
node  number 
generation  increment 
DOF ( I )  force  ( di spl . ) 
increment  factor  for  I-force 
( di spl . ) 


Terminate  with  a  blank  card.  Boundary  data  generation  algorithm  is  iden 
tical  to  the  node  geometry  generation  (see  Note  2). 


4.3.7  Material  and  Element  P roperty  Card 

4. 3. 7.1  Macro  Control  Card 

The  macro  instruction  MATE { A4 )  AT  THE  FIRST  CARD  causes  material  data 
input  inode. 

4. 3. 7. 2  Material  Set  Heading  Card  (20A4) 

A1 pha-numeric  heading  indicating  the  material,  element  type  and  so  on 
is  given.  One  card  must  be  provided  for  each  material  set. 
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4. 3. 7. 3  Identification  Cards 

NAMELIST  input  with  the  name  MAUD  at  commencement. 


Vari  able 
MAT  NO 
NGAUSS 


NSTRES 


Desc  ri ption 
Material  set  number 

Humber  of  Gauss  integartion  points  (1  to  4). 

Need  not  to  be  provided  if  recommended  default  value 
2  is  to  be  used. 

Number  of  stress  points  (1  to  4).  Need  not  to  be  pro¬ 
vided  if  recommended  defeult  value  of  1  is  used. 


LTYPE  Element  type  for  this  material  set; 

9  for  plane  strain  consolidation, 

10  for  one-dimensional  dynamic  analysis, 

11  for  dynamic  analysis  of  plane  strain  case. 


Example  of  identification  card; 


bATIDb  MATN0=1 ,NGAUSS=1 ,LTYPE=10 


4. 3. 7. 4  Material  Property  Card 

NAMELIST  input  with  the  name  MATRL (A7 )  at  commencement. 

(a)  Material  properties  required  for  quasi-stati c  consol i dati on 
Variable 
El 
E2 

XNU 1 
XNU2 
G2 
RHO 


Descri  ption 

Modulus  of  elasticity  1-axis. 

2-axis. 

Poisson's  ratio  in  1-  and  2-axis. 

1-  and  3-axis. 


PHI 


Shear  modulus  (to  be  computed  if  not  given). 
Mixture  mass  density. 

Reciprocal  of  fluid  compressibility. 


K 1 


Permeability  in  1  - di recti  on. 

2-di recti  on . 
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K2 

Cl 


FI ui d  mass  densi ty . 


(b)  Material 
Variable 
El 
E2 

XNU1 

XNU2 

G2 

RHO 

PHI 

Kl 

K2 

Cl 

ALPA1 

ALPA2 

SC 

ST 


properties  required  for  dynamic  analysis 
Oescri ption 

Modulus  of  elasticity  1-axis 

2-axi s 

Poisson's  ratio  1-  and  2-axis 

1-  and  3-axi s 

Shear  modulus  (to  be  computed  if  not  given) 
Mixture  mass  density  (?) 

Porosi  ty 

Permeabi 1 i ty  in  1- direction 

2-  di recti  on 

Fluid  mass  density  C?2) 

Fluid  compressibi 1 i ty 
Solid  compressibility 
Stiffner  damping  coefficient  for  solid 
Mass  damping  coefficient  for  solid 


4.3.8  Macro  Commands  for  Solution  Phase 

For  solution  process,  the  following  cards  are  needed. 

(a)  For  Quasi-static  Consolidation  Analysi s 

Col.  1-4  Col.  6-9  Description 

TANG  BAN  Symmetric  tangent  matrix 


FORM 

BFO 

Form  body  force  vector 

__  ★  ★ 

CONS 

SAVET 

Perform  consolidation  analysis 

STOP 

Terminate  the  run 

k  If  body  force  is  not  considered,  do  not  include  this  card.  However, 
if  this  option  is  used,  1-  and  2-components  of  the  gravi tational 
constant  with  correct  signs  must  be  provided  in  (2F10.0)  format  im¬ 
mediately  after  'FORM  BFO'  card. 

**  If  this  option  is  used,  the  output  will  be  saved  on  NTAPE  (NTAPE=11 
for  displacement,  12  for  stresses). 

***  This  card  must  be  followed  by  the  following  cards. 

1st  card  (3F10.0) 


Col umns 
1-10 
11-20 
21-30 


Vari  able 
Uo 
Vo 
Po 


Description 

Initial  horizontal  displacement 
Initial  vertical  displacement 
Initial  pore  pressure 


2nd  Card  Set  (F_1_0.0,4I_5 ) 

Repeat  as  many  NSTEPS  as  defined  in  C.4.2. 

1-10  (TIMES(I ) ,1=1, NSTEPS )  Time  at  the  end  of  a  step 
11-15  (MDT(I),  1  =  1, NSTEPS)  Number  of  time  increments/step 

16-20  ( NOP T ( I ), 1  =  1 .NSTEPS )  Time  domain  integration  factor 

21-25  (NPRINT( I ) , I =1 , NSTEPS )  Printout  parameter 
26-30  (KPRS(I),  1=1, NSTEPS)  =0/1  do  not  print/print  stress 

where  NOPT  =  1  for  d=0.5 
=  2  for  o(,=0.67 
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=  3  for  ot=0 .871 


=  4 

for  Or  1.0 

=  5 

for  oH  +  l/At 

- 1 / 1 n ( 1 +  At  > 

(b)  Dynamic  Response  Analysis 

Col.  1-4 

Col. 6-9 

Descri ption 

TANG 

BAN 

Symmetric  tangent  matrix 

FORM* 

BFO 

Form  body  force  vector 

_  ★★★ 

DYNM 

★  ★ 

SAVET 

Perform  dynamic  response  analysis 

STOP 

Terminate  the  run 

*  Same  as  in  Part  a 
**  Same  as  in  Part  a 

***  This  card  must  be  followed  by  the  following  cards 
1st  Card  (4F1Q.0) 


Col. 

Variable 

Description 

1-10 

uso 

I ni tial 

solid  displacement 

11-20 

vso 

Initial 

solid  velocity 

21-30 

UFO 

Initial 

fluid  relative  displacement 

31-40 

VFO 

I  ni ti al 

fluid  relative  velocity 

2nd  Card  (I5.4F10.0) 

The  dynamic  load  is  input  by  the  coefficients  of  loading  function 
defined  as 

P(t)  =  AZERO  +  RZERO*COS ( 0MEGA*T )  +  CZERO*SIN(OMEGA*T) 

Col.  Variable 


1-5 


LOAD 


Oescri ption 
Loading  criterion-1" 


6-15 


OMEGA 


16-25  AZERO 
26-35  BZERO 
36-45  CZERO 


Loading  Frequency  (radius/sec) 

Coef fici ent 

Coefficient 

Coefficient 


+  LOAD  =  0  ;  Loading  by  traction  only 

=  1  ;  Loading  by  support  accel leration  only 
=  2  ;  Loading  by  traction  and  acceleration 

If  LOAD  =  0,  skip  the  following  two  cards  and  go  to  5th  card. 

3rd  Card  (415) 

Col.  Variable  Description 

1-5  LACEL  Acceleration  with  time  ++=0/l 

6-  5+NDF*5  IACEL  =0/1  for  which  degree  of  freedom 

has  acceleration 

++  LACEL=0  Use  the  loading  function  in  Eq.(B-l) 

LACEL=1  Use  the  Accel egram  recorded  data 

If  LACEL  =0,  skip  the  following  card  and  go  to  5th  card. 

4th  card  ( 215) 

Col.  Variable  Description 

1-5  NTAPE  The  tape  No.  on  which  accelegram  data 

is  given 

6-10  NDTACL  Number  of  acceleration  records 

Then,  the  following  data  are  read  from  the  tape; 

(TREC0RD( I ) ,ACEL ( I ) , 1=1 .NDTACL ) 
where  TRECORD  is  the  time  and  ACEL  is  the  acceleration  value. 
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5th  Card  (4F10. 0,315) 

Repeat  as  many  as  NSTEPS  given 

Col.  Variable 

I- 10  (Times ( 1 ) , 1=1 .NSTEPS) 

II- 20  (0( I ) ,1=1 .NSTEPS) 

21-30  (p(I), 1=1, NSTEPS) 

31-40  (YU),  1  =  1,  NSTEPS) 

41-45  (NDT(I), 1=1, NSTEPS) 

46-50  (NPRINT ( I ) , 1+1 , NSTEPS) 
51-55  (KPRS(I), 1=1, NSTEPS) 


in  Control  Data  Card. 

Descri  ption 

The  time  at  the  end  of  a  step 
Wilson' s  coefficient 
Newmark's  coefficient 
Newmark's  coefficient 
Number  of  increments/steps 
=J;  Print  Jth  increment 
=0/1  Print/do  not  print  stresses 
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Appendix  A 
PROGRAM  LISTING 


The  complete  listing  of  the  computer  program  is  given  below. 

_ PROGRAMMED  BY  PROF.  J.  K.  LEE,  B.L .  ABOUSTIT 

DEPT  OF  ENGR.  MECH.OHIO  STATE  UNIVERSITY,  COLUMBUS,  OHIO. 

_ MAIN  PROGRAM  :  SET  PROGRAM  CAPACITY. 


....PROBLEM  SIZE  IS  CONTROLLED  BY  THE  DIMENSION  IN  BLANK  COMMON  AN 
THE  VALUE  OF  MAX. 

_ ALL  ARRAYS  RESIDE  IN  CENTRAL  MEMORY. 

_ DOUBLE  PRECISION  IS  ASSUMED.  FOR  SINGLE  PRECISION  SET  IPR  =  1 

AND  SUPPRESS  ALL  IMPLICIT  REAL  STATEMENTS. 

CALLED  BY  - 
CALLS  -  ACNTRL 

....  FILE  ASSIGNMENT  : 

FILE  NO.  10  =  MESH  DATA  IF  NODE  SAVET,  ETC. 

FILE  NO.  11  =  SOLUTION  (DISP)  IF  DISP  SAVET. 

FILE  NO.  12  =  STRESSES  IF  STRE  SAVET 
FILE  NO.  13  =  SCRATCH  FILE  IF  NTEMP  .GE.  2. 

FILE  NO.  14  =  SCRATCH  FILE  FOR  STRAINS 

FILE  NO.  15  =  SCRATCH  FILE  FOR  ELEMENT  STIFFNESS 

FILE  NO.  16  =  SCRATCH  FILE  FOR  GLOBAL  STIFFNESS 

FILE  NO.  17  =  SCRATCH  FILE  FOR  GLOBAL  FORCE  VECTOR 

FILES  10  TO  12  USE  FORMATTED  I/O  (  5  FIELD  INTEGERS  AND  12  FIELD 
REALS  ). 

COMMON  M( 280000) 

COMMON  /SIZE/  MAX 
MAX=280000 
IPR=2 
C 

CALL  ACNTRL  (IPR) 

C 

STOP 

END 
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BLOCK  DATA 


IMPLICIT  REAL*8(A-H ,0-Z ) 

COMMON  /CDATA/  0 .HEAD ( 20 ) .NUMNP .NUMEL .NUMMAT.NEN ,NEQ .MEQ.NTAPE ,NPN 
COMMON  /GQDRT /  GPT(4,4) ,WT{4 ,4)  ,SHP (4 ,9 )  ,LX (4 ) ,LE (4) 

COMMON  /LABL/  PD (9 ) ,A{9 ) ,BC (2 ) ,01 (6 ) ,CD( 3 ) ,TE ( 3 ) ,FD{ 3 ) ,FOL (3 ) ,FON ( 
14) 

COMMON  /ANGLE/  ANG (4) , DC ( 3 ) , ST ( 2 ) , IEPS ,L INE ,TWR 

COMMON  /MACRO/  WD (15) .NWD.NV .NC.NRD 

COMMON  /ML  I  ST /  IEL .NTEMP .NGAUSS .MSTRES ,DUMY (23) 

COMMON  /ELABL/  ITYPE(6,10) 

DATA  F0L/3H ( 1 5  ,2H  ,4HI5)  /,F0N/3H( I5.2H  .4HF12 . ,2H4 )/ 

DATA  ANG/57 .29578,6.2831853 ,32. 174, 9. 81/, ST/3HESS.3HAIN/ 

DATA  A/2H,1,2H,2,2H,3,2H,4,2H,5,2H,6,2H,7,2H,8,2H,9/,CD/4H  C00.4HR 
1DIN.4HATES/ 

DATA  TE/4M  TEM.4HPERA ,4HTURE/,FD/4H  B0U.4HND  V.4HALUE/ 

DATA  CC/4H  B.C.2H.  /.DI/4H  DIS.2HPL.4H  VEL,2H0C,4H  ACC.2HEL/ 

DATA  0/1H0/ 

DATA  PD/6H ( 14 ,2X ,2H  ,2HI1,2H  .8HD12.5  ,2H  ,8HD12.5,2X,2H  ,6H 

1D12.5)/ 

DATA  GPT/4*0.D0, -.577350269189600, .577350269 189600,2*0. DO, -.774596 
16692415D0.0. DO,. 774596669241500,0. DO, - . 861 1 36 3 1 1 5941  DO , - . 339981043 
25849D0, .3399810435849D0, .8611363115941D0/ 

DATA  WT /2 .DO ,3*0. DO ,2*1 .DO ,2*0.00 , . 5555555555556D0 , .8888888888889D 
10 , . 555555 5555556D0, 0. DO, .3478548451 3 75D0 ,  .652 1451 548625D0 , .6521451 
2548625D0 , .347854845 13 75DO/ 

DATA  WD/4HN0DE ,4HELEM,4HB0UN .4HMATE  ,4HBVAL ,4HTEMP ,4HUTAN ,4HTANG ,4H 
1F0RM,4HS0LV ,4HDISP .4HSTRE ,4HC0NS ,4HDYNM,4HST0P/,NWD/15/,NV/1/,NC/1 
2/.NRD/26/ 

DATA  LX/-1 ,1 ,1 ,-l/,LE/-l . -1,1,1/ 

DATA  NTEMP/1 /,NGAUSS/2/,NSTRES/l/,DUMY/23*0.0D0/ 

DATA  ITYPE/4HPLAN ,4HE  ST ,4HRESS,3*4H  ,4HPLAN,4HE  ST ,4HRAIN ,3*4H 
1  ,4HAXIS,4HYMME,4HTRIC,4H  STR.4HESS  ,4H  ,4H1-D  ,4HCART ,4HESI 

2A.4HN  TR.4HANSP ,4H0RT  ,4H2-D  ,4HCART ,4HESIA,4HN  TR .4HANSP .4H0RT  ,4 
3H1-D  .4HAXIS ,4HYM.  ,4HTRAN ,4HSP0R ,4HT  ,4H2-D  .4HAXIS ,4HYM.  ,4HTR 
4AN .4HSP0R ,4HT  ,4H2-D  , 4HADVE , 4HCT 1 0 , 4HN-D 1 , 4HFUSS , 4H I ON  ,4HP.ST, 
54HRAIN ,4H-C0N ,4HS0L I ,4HDATI ,4H0N  ,6*4H  / 

END 
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SUBROUTINE  ACNTRL  (IPR) 

C  - - - 

IMPLICIT  REAL*8(A-H ,0-Z) 

REAL *8  NP 

....MAIN  CONTROL  ROUTINE 

_ READ  HEADING  AND  THE  CONTROL  DATA. 

....SET  STORAGE  POINTERS 

CALLED  BY  -  MAIN 

CALLS  -  PCOMP .ZERO ,GENVEC .CPTABL .BCCODE ,PR0FIL , 

setmpt,eldata,elmlib,pload,formel,actcol 

UACTCL.PRTDIS, CONSOL 

COMMON  /CDATA/  0  .HEAD (20 )  .NUMNP  .NUMEL  .NUMMAT ,NEN  ,NEQ  ,MEQ  .NTAPE  ,NPN 
COMMON  /L ABL /  PD (9 )  ,A( 9 )  ,BC (2 ) ,DI (6 ) ,CD ( 3 ) , TE ( 3 ) ,FD (3 ) ,FOL ( 3 ) ,FON ( 
14) 

COMMON  /PRLOD/  PROP 

COMMON  /ANGLE/  ANG (4 ) ,DC ( 3 ) , ST( 2  ) , IEPS ,L I NE , TWR 
COMMON  /TIMECK/  CPTIME (20 ) .COMPUT (20 ) ,NCHK 
COMMON  /MACRO/  WD (15) ,NWD ,NV ,NC ,NRD 

NAMELIST  /CONTRL/  NUMNP .NUMEL .NUMMAT ,NEN .NDM.NDF ,NPN ,N STEPS 
LOGICAL  PCOMP .ERR.PRT .UNSYM.AFAC .BACK ,AFL ,BFL ,CFL , TRD ,TWR ,TLD ,BND 
COMMON  M ( 1 ) 

DIMENSION  TI TL (20 ) ,  SAVEFL(lO),  NFILE(IO) 

DATA  START/4HLAMP/.C0MNT/4HC. . . / , TL /4H  TL0/.BF/4H  BF0/.NP/4H  NOP/, 
1R0/4H  REA/.WR/4H  SAV/.PR/4H  PRI/.BN/4H  BAN/ 

_ FILE  NO.  10  :  MESH  DATA 

....FILE  NO.  11  :  SOLUTION  (  IN  SUB  PRTDIS  ) 

_ FILE  NO.  12  :  P. STRESSES  (  IN  SUB  PVALUE  ) 

ERR=. FALSE. 

PR0P=1 .ODO 

C  ....READ  AND  WRITE  INPUT  DATA  AS  GIVEN 

LIME  =0 

WRITE  (6,38) 

1  READ  ( 5 ,34  ,END=2 )  TI TL 
LINE=LINE+1 

WRITE  (6,39)  LINE  , T I TL 
GO  TO  1 

2  REWIND  5 
NSAVE=0 

3  IF(ERR)  STOP 
READ  (5,34)  TITL 
WRITE  (6,40)  TITL 

IF (PCOMP (TI TL ( 1 ) ,COMNT ) )  GO  TO  3 
DO  4  I M= 1 ,NWD 

IF (PCOMP (TITL ( 1 ) ,W0 ( I M ) ) )  GO  TO  6 

4  CONTINUE 

IF ( . NOT. PCOMP (TI TL ( 1 ), START)  )  GO  TO  31 
DO  5  1=1,20 

5  HEAD ( I ) =TITL ( I ) 

IF(NCHK.NE.O)  CALL  CPTABL  ( TI TL ( 1 ) ,  1  ) 

CALL  ZERO  (CPTIME, 20) 


o  o  <~>  o  o  o  o  o  o 
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NCHK-0 

c 

READ  ( 5  ,CONTRL  ) 

WRITE  (6,37)  HEAD ,MUMNP ,NUMEL ,NUfWAT ,NEN ,NDM,NDF 


C  NUMNP  NUMBER  OF  NODAL  POINTS 

NUMEL  NUMBER  OF  ELEMENTS 

NEN  MAX.  HUMBER  OF  NODES  PER  ELEMENT 

NDM  NUMBER  OF  DIMENSION 

NDF  NO.  OF  DEGREE  OF  FREEDOM  PER  NODE 


MEQ=NUMNP*NDF 
NST=NDF*NEN 
NEN1=NEN+1 
MEQR=MEQ*IPR 
AFAC= .TRUE . 

BACK=.TRUE. 

TLD=. FALSE. 

....SET  STORAGE  POINTERS  FOR  ALL  ARRAYS  EXCEPT  FOR  THE  GLOBAL  STIF 
_ M'S  ARE  FOR  INTEGER  ARRAYS  AND  N'S  ARE  FOR  REAL  ARRAYS. 

NO-1 

NT=NO+MUMNP*NDM*IPR 
MO-NT +NUMNP*IPR 
Ml=MO+NENl*NUMEL 
M2=M1+NUMMAT*4 
M3=M2+NDF*NUMNP 
M4-M3+NST 
M5-M4+NEN 
LAST-M5 

CALL  SETMPT  (N1 ,M5,MEQ,LAST,TITL(1 ) ) 

N2=N1+NDM*NEN*IPR 
N3=N2+NRD*IPR 
N4=M3+MST*MST*IPR 
N5=N4+NST*IPR 
N6=N5+NST*IPR 
N7=N6+NEM*IPR 
N8-N7+MEQR 
N9-N8+MEQR 
N1O-N9+ME0R 

CALL  ZERO  (M(N7),MEQ*3) 

GO  TO  3 

C  CHECK  MACRO  COMMAND  CARD  AND  TRANSFER  TO  APPROPRIATE  PROCEDURE 
6  TRD=PCOMP(TITL(2) ,RD) 

TWR-PCOMP (TITL (2 ) ,WR ) 

PRT-.TRUE . 

IF( IM.LE.7)  NTAPE-10 

IF( IM.EQ.ll . OR .IM.EQ.13)  NTAPE=11 

IF ( IM.EQ. 12 )  NTAPE-12 

I F ( PCOMP ( T I TL (2 ) ,NP  ) )  PRT-. FALSE. 

IF(.NOT.TRD)  GO  TO  8 
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IF (NCHK.NE .0)  GO  TO  7 

READ  ( NTAPE  ,36  )  TITL.  .NODES  ,NELS  ,NUMS  ,NENS  ,NDMS,NDFS 
IF (NODES. EQ.NUMNP. AND. NELS.EQ.NUMEL. AND. NEWS. EO.NEN. AND. NDMS.EQ.ND 
1M.AND.NDFS.EQ.NDF)  GO  TO  7 

WRITE  (6,46)  NTAPE, NODES, NELS.NUMS, WENS, NDMS,NDFS 
STOP 

7  READ  (NTAPE, 34)  TITL(l) 

IF (PCOMP (TITL ( 1 ) ,WD ( IM) ) )  GO  TO  9 
WRITE  (6,47)  TITL, WD(IM), NTAPE 
STOP 

8  IF ( .MOT.TWR)  GO  TO  9 

IF(NCHK.EQ.O)  WRITE  (NTAPE, 36)  HEAD  .N'JMNP  ,NUMEL  .NUMMAT  ,NEN  ,NDM,NDF 
WRITE  (NTAPE, 34)  TITL (1 ) 

NSAVE=NSAVE+1 
SAVEFL ( NSAVE )=TI TL ( 1 ) 

NF1LE (NSAVE )=NTAPE 


9  CALL  SCL0K1 

GO  TO  (10,11,12,13,14,15,16,17,22,24,27,28,29,30,33),  IM 

C  . READ  AND  /OR  GENERATE  NODAL  GEOMETRY  DATA _ 

C  TITL ( 1 )  -  NODE 

10  CONTINUE 

CALL  GENVEC  ( NDM ,M( NO ) ,CD ,PRT ,ERR ,TRD ,TWR ) 

GO  TO  32 

C  . READ  AND/OR  GENERATE  ELEMENT  DATA . 

C  TITL ( 1 )  =  ELEM 

11  CONTINUE 

CALL  ELDATA  (M(MO) ,M(M4 ) ,NEN1 ,NDF ,PRT ,ERR ,TRD ,TWR ,MD IF ) 

GO  TO  32 

C  READ  AND/OR  GENERATE  BOUNDARY  CONSTRAINT  CODES. 

C  TITL ( 1 )  =  BOUN 

12  CALL  BCCOOE  (M(M2) ,M(M3) ,NDF,NUMNP,PRT, TRD,TWR, HEAD, NTAPE ) 
GO  TO  32 

C 

C  _ READ  THE  MATERIAL  DATA. 


C 

13  WRITE  (6,41)  HEAD 
N0W=0 

CALL  MATLIN  (M(N10 ) ,NADD ,M { Ml ) .NUMMAT ,NRD  ,M(N2 ) ) 

LAST=N10 

CALL  SETMPT  (Nil ,N10 ,NADD*IPR, LAST  ,TI TL ( 1 ) ) 

GO  TO  3 
C 

C  _  GET  NODAL  30UNDARY  VALUES  (  NATURAL  AND/OR  ESSENTIAL  ). 

C 

14  CALL  GENVEC  (NDF ,M(N8)  ,FD ,PRT,ERR,TRD  ,TWR) 

GO  TO  3 

C 

C 

15  CALL  GENVEC  ( 1 ,M( NT ) , TE ,PRT , ERR.TRD ,TWR ) 

GO  TO  3 

C 

C  ....  FORM  TANGENT  MATRIX 


O  O  £->  O  O  O 


16  UNSYM=.TRUE. 

17  IF( IM.EQ.8)  UNSYM= .FALSE . 

CFL=UNSYM 

AFL=.TRUE. 

BFL=.TRUE. 

AFAC=.TRUE. 

BND=PCOMP (TITL (2 ) ,BN ) 

IF ( .N0T.8ND  )  GO  TO  18 
I BAND= ( MD I F+l )*NDF 
WRITE  (6,50)  I  BAND 
NCOEF=MEQ*IBAND 
GO  TO  19 

18  CALL  PROFIL  (M(M5 ) ,M(M2 ) ,M(MO) ,NDF,NEN1 ,NCOEF ,TRD ,TWR ) 

19  CALL  SETMPT  (N12 ,N11 , NCOEF * I P R , LAST  ,TITL ( 1 ) ) 

MA=N1 1 

N 1 3=N 1 2 

CALL  ZERO  (M(NA).NCOEF) 

IF( .NOT.UNSYM)  GO  TO  20 

CALL  SETMPT  (N13 ,N12 ,NCOEF*IPR .LAST , TITL ( 1 ) ) 

NC =N 1 2 

CALL  ZERO  (M(NC).NCOEF) 

20  WRITE  (6,44)  NO, NT ,M0,M1 ,M2 ,M3 ,M4,M5,N1  ,N2 ,N3 ,N4,N5 ,N6 ,N7 ,N8,N9 ,N1 
10 ,NA,NC  ,N13 

WRITE  (6,45)  MEQ ,NEQ ,NCOEF .LAST ,NRD 
ISVJ=2 

21  N0W=N0W+1 
PRT=PCOMP (TITL (3) ,PR ) 

CALL  FORMEL  (M(NO) ,M(N1 ) ,M(N10 ) ,M(N2) ,M(N3) ,M(N4) ,M(N5 ) ,M(N6 ) ,M(N7 
1 ) ,M(N8) ,M(NT) ,M(MO) ,M(M1 ) ,M(M2) ,M(M3 ) ,M(M4) ,M(M5) ,NDF,NDM,NEN1 ,NST 
2  ,ISW,M(NA ) ,M(NC ) ,M(N9) ,AFL ,BFL ,CFL ,PRT,TLD ,NRD ,NOW, IBAND ,BND ) 

IF ( ISW.NE . 5 )  GO  TO  32 
TLO=. FALSE. 

CALL  ZERO  (M(N7),MEQ) 

GO  TO  32 

-  FORM  CONSISTANT  LOADING  DUE  TO  TEMPERATURE  OR/AND  BODY  FORCE. 

22  ISW=0 

I F ( .NOT.PCOMP (TITL (2) ,TL ) )  GO  TO  23 
ISW=3 

TLD=.TRUE . 

23  IF(PCOMP (TITL (2 ) ,BF) )  ISW=4 
IF( ISW.EQ.O)  GO  TO  3 

IF ( ISW.EQ.4)  READ  (5,35)  ( DC ( I ) , 1=1 ,NDM) 

- DC ( I )  =  I  COMPONENT  OF  THE  GRAVITATIONAL  CONSTANT  G. 

AFL=. FALSE. 

CFL=. FALSE. 

BFL= .TRUE . 

PRT=. FALSE. 
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GO  TO  21 
_  SOLVE  THE  EQUATIONS. 

24  CALL  PLOAD  (M(M2 ) ,M(N8 ) ,M(N7 ) ,M(N9 ) ,M(NT ) ,MEQ , 1 .DO ,0,BND ) 

IF ( UNSYM )  GO  TO  25 

IF(BND)  CALL  BANSOL  (M(MA ) ,M(N7 ) ,M(N7 ) ,MEO , I  BAND ,AFAC .BACK ) 

IF ( .NOT.BND  )  CALL  ACTCOL  (M( NA ) ,M(N7 ) ,M(M5 ) ,MEQ , AFAC .BACK ) 

GO  TO  26 
CONTINUE 
AFAC=. FALSE. 

IF (PCOMP (TI TL ( 2 ) ,TE ( 1 ) )  )  CALL  PLOAD  ( M(M2 ) ,M (N8 ) ,M( N7 ) ,M( N9 ) ,M(NT ) 
1  ,MEQ , 1 .ODO  ,1  ,BND  ) 

CALL  ZERO  ( M ( N 9 ) ,MEQ) 

GO  TO  32 

....  PRINT  THE  SOLUTION  VECTOR 

CALL  PRTDIS  (M(M2) ,M(NO) ,M(N7) ,M(N8) ,NDM,NDF,TRD .TWR.PRT, 0,0,0. DO, 
1BND ) 

GO  TO  3 

....  COMPUTE  AND  PRINT  GRADIENT  OF  SOLUTION  (  OR  STRESS  ) 

AT  GAUSS  PTS . 


ISW=5 

AFL=. FALSE. 

BFL=. FALSE. 

CFL=. FALSE. 

LINE=0 
GO  TO  21 

_  CONSOLIDATION  ANALYSIS  _  TITL(l)  =  CONS  _ 

CONTINUE 

M ( N 1 3 )  =  TIMES! 1 ) 

M(M6)  -  NDT(l),  M(M7)  =  NPRINT(l),  M(M8)  =  NOPT(l) 

M6=N13+NSTEPS*IPR 

M7=M6+NSTEPS 

M8=M7+NSTEPS 

CALL  SETMPT  (N14 ,M8 ,NSTEPS .LAST  ,TI TL ( 1 ) ) 

WRITE  (6,48)  N13 ,M6 ,M7 ,M8,N14,LAST 

CALL  CONSOL  (M(NO ) ,M(NT ) ,M(MO ) ,M(M1 ) ,M(M2 ) ,M(M3 ) ,M(M4 ) ,M(M5 ) ,M(N1 ) 
l,fl(N2) ,M(N3) ,M(N4) ,M(N5) ,M(N6) ,M(N7) ,M(N8) ,M(N9) ,M(N10) ,M(NA) ,M(N1 
23  )  ,M(M6 ) ,M(M7),M(M8) ,NEN ,NEN1 ,NPN ,NDF ,NDM,MEQ , IBAND.NUMEL ,NUMNP ,NS 
3TEPS ,NST ,NRD ,BND ) 


GO  TO  32 
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_  DYNAMIC  CONSOLIDATION  ANALYSIS 

CONTINUE 

M { N 1 3 )  =  TIMES! 1 ) ,  M(M6)  =  NDT(l),  M(M7)=  NPRINT(l) 
M(N14 )  =  THETA! 1 ) ,  M(N15)  =  GAMA ( 1 ) ,  M(N16)  =  BETA! 1 ) 

M  ( N 1 7  )  =  VO(NDF , 1 ) ,  M(N18)  =  AO(NDF,l) 

M( N19 )  =  U ( NDF , 1 ) ,  M(N20)  =  V (NOF.l ) ,M(N21 )  =  A(NDF,1 ) 


M6=N13+NSTEPS*IPR 

M7=M6+NSTEPS 

LAST=M7 

CALL  SETMPT  (N14 ,M7 ,NSTEPS  .LAST .TITL ( 1 ) ) 

N15=N14+NSTEPS*IPR 

N16=N15+NSTEPS*IPR 

N17=N16+NSTEPS*IPR 

N18=N17+NDF*NUMNP*IPR 

N19=N18+NDF*MUMNP*IPR 

N20=N19+NDF*NUMNP*IPR 

N21=N20+NDF*NUMNP*IPR 

LAST=N21 

CALL  SETMPT  (N22.N21 , NDF *NUMNP*IPR, LAST, TITL(1 ) ) 

PRINT  49,  N13,M6,M7,N14,N15,N16,N17,N18,N19,N20,N21,N22,LAST 

CALL  DYNAMC  (M( NO ) ,M(NT ) ,M(MO ) ,M!M1 ) ,M(M2 ) ,M(M3 ) ,M(M4) ,M(M5 ) ,M(N1 ) 
1  ,M(M2 ) ,M(N3 ) ,M(N4) ,M(N5 ) ,M(N6) ,M(N7 ) ,M(N8 ) ,M(N9 ) ,M(N10 ) ,M(NA) ,M(N1 
23 )  ,M(M6 ) ,M(M7 ) ,M(N14 ) ,M(N15 ) ,M(N16 ) ,M(N17 ) ,M(N18 ) ,M(N19 ) ,M(N20) ,M( 
3N21 )  ,NEN ,NEN1 ,NDF,NDM,MEQ , IBAND .NUMEL ,NUMNP ,NSTEPS,NST ,NRD,BND ) 

GO  TO  32 

WRITE  (6,42)  TITL 
STOP 

CALL  CPTABL  (TI  ' 1 )  ,0 ) 

GO  TO  3 
CONTINUE 

CALL  CPTABL  (TITL ( I ) ,1 ) 

IF ( NSAVE . NE .0 )  WRITE  (6,43)  (N  ,NF I LE ( N  )  .SAVEFL (N ) ,N  =  1 ,NSAVE ) 

RETURN 


FORMAT  ( 20A4 ) 

FORMAT  (3FI0.0) 

FORMAT  (20A4 ,/ ,1015) 

FORMAT  (1H1,5X,20A4,//,5X, 'NUMBER  OF  NODAL  POINTS 

1  =  '  ,15,/, 5X, 'NUMBER  OF  ELEMENTS 

2  =  '  ,15,/, 5X, ‘NUMBER  OF  MATERIAL  SETS 

315, /,5X, 'MAXIMUM  NUMBER  OF  NODES  PER  ELEMENT 

4X, 'DIMENSION  OF  COORDINATE  SPACE 

5REE  OF  FREEDOM  PER  NODE  =  ' , I 

FORMAT  ( 1111  ,//10X  , ' _ LIST  OF  INPUT  AS  GIVEN _ 

FORMAT  (5X , ' CARD  NO.  =' ,14, ' : ' ,20A4) 

FORMAT  (//5X.20A4) 


=  M5,/,5 
=  ' ,I5,/,5X,'DEG 
,15,//) 

..'//) 
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41  FORMAT  (1H1.5X.20A4) 

42  FORMAT  (///5X,‘ _ FATAL  ERROR _ THE  MACRO  COMMAND '//5X ,20A4//5X , 

1  *  IS  MOT  ALLOWED  IN  THIS  PROGRAM.') 

43  FORMAT  (//// lQX.'D  A  T  A  FILE  * //10X , 'FILE  NO. ' ,5X, 'MACRO  7/10 
1 (5X ,215 ,5X ,A4/) ) 

44  FORMAT  ( ///5X , ‘ ARRAY  POINTERS  :  NO ,NT,M0-M5 .N1-N10 ,NA ,NC ,N13 7/5X , 
12(12110//) ) 

45  FORMAT  ( ///5X NUMBER  OF  TOTAL  EQUATIONS  (MEQ  )  =', I8/5X NUMBER  OF 
1  ACTIVE  EQUATIONS(NFQ)  =', I8/5X NUMBER  OF  COEFFICIENTS  IN  (A) 

2' ,I8/5X, 'REQUIRED  SIZE  OF  (M)  = ’ , I8/5X , ’ MRD  (  SET  IN  3L 

30CK  DATA  )  =  * , 18) 

46  FORMAT  (///5X,1 _ FATAL  ERROR _ CONTROL  DATA  FROM  TAPE  NO. ',15, 

1’  ARE  : ' ,//5X , ' NUMBER  OF  NODAL  POINTS  =  1 

2, 15, /,5X, 'NUMBER  OF  ELEMENTS  =  ',15,/, 

35X, 'NUMBER  OF  MATERIAL  SETS  =  ' , 1 5 , / ,5X , ' MA 

4X1 MUM  NUMBER  OF  NODES  PER  ELEMENT  =  ' , 15 ,/,5X , ' D IMENSIO 

5N  OF  COORDINATE  SPACE  =  ' ,15,/, 5X , 'DEGREE  OF  FREE 

600M  PER  NODE  =  ' ,15,/, 5X,//,5X, 'WHICH  DO  NOT  M 

7ATCH  WITH  THE  CONTROL  DATA  SPEC IFIEDD  FOR  THIS  RUN .  EXC 

8UTI0N  TERMINATED .  ') 

47  FORMAT  ( ///5X , ' . FATAL  ERROR..  :  MACRO  MISMATCH . CORRECT 

1  MACRO  IS  ' ,//5X,20A4,//5X, 'INSTEAD  OF  *',A4,'*  TO  RETRIEVE 

2  DATA  FROM  THE  TAPE  NUMBER ',15,' _ EXCUTION  TERMINATED _ ') 

48  FORMAT  (///5X , 'ADDITIONAL  MEMORY  POINTERS  FOR  THE  CONSOLIDATION  AN 

1ALYSIS ' //5X , 'TIMES(1 )  =  M(N13)  : ' , I8/5X , ' NDT ( 1 )  =  M(M6) 

2  : ' , I8/5X , ' NOPT ( 1 )  =  M(M7)  : ' , I8/5X , ' NP  R I  NT ( 1 )  =M(M8)  : 

3 ' , I8/5X , ‘ TOTAL  RESERVED  M(N14)  : '  ,  18 ,5X , ' LAST  =',I8) 

49  FORMAT  (///5X , 'ADDITIONAL  MEMORY  POINTERS  FOR  THE  DYNAMIC  ANA 

1  LYSIS'/ /5X , 'TIMES ( 1 )  =  M(N13)  : '  ,  I8/5X , 'NDT ( 1 )  =  M(M6) 

2  : ' ,I8/5X, 'NPRINT(l)  =M(M7)  : ' ,I8/5X , ' THETA ( 1 )  =  M(N14)  : 

3 ' ,  1 8 /5X  , '  GAMA  ( 1 )  =  M  ( N 1 5  )  :  '  ,I8/5X, 'BETA(l)  =  *1  ( N 16  )  :', 

4I8/5X  , ' VO ( 1 )  =  M ( N 1 7 )  : '  ,  I8/5X  , ' AO ( 1  )  =M(N18)  :  '  ,  18 

5/5X , ' U ( 1 )  =  T 1  ( N 1 9  )  : ' , I8/5X , ' V ( 1 )  =  M(N20)  : * , 13/5 

6X  , ' A( 1 )  =  M(M21)  , I8/5X , 'TOTAL  RESERVED  M(N22)  :' ,18, 5X  , 

7 ' LAST  =' ,18) 

50  FORMAT  ( //5X , ' BAND  SOLVER  IS  USED . IBAND  =’  ,18) 

END 


SUBROUTINE  ACTCOL  (A,B ,JD IAG ,NEQ,AFAC .BACK ) 


IMPLICIT  REAL  *8 ( A-H ,0-Z  ) 

LOGICAL  AF AC, BACK 
COMMON  /ENGYS/  AENGY 
DIMENSION  A(l),  B(l),  JDIAG(l) 

....  ACTIVE  COLUMN  PROFILE  SYMMETRIC  EQUATION  SOLVER 

CALLED  BY  -  ACNTRL 
CALLS  -  DOT 

....  FACTOR  A  TO  UT*U*U,  REDUCE  B 

AENGY=O.ODO 

JR=0 

DO  6  J  =  1  ,NEQ 
JD=JDIAG(J ) 

JH=JD-JR 

I S  = J -JH+2 

IF ( JH-2 )  6,3,1 

IF ( .NOT.  AF'.C )  GO  TO  5 

I E  =J  - 1 

K=JR+2 

ID=JD IAG ( I S - 1 ) 

_  REDUCE  ALL  EQUATIONS  EXCEPT  DIAGONAL 

DO  2  I=IS, IE 
IR=ID 

I  D=JD  I  AG  ( I ) 

IU=MINO( ID-IR-1 , I -IS+1 ) 

IF ( IH.GT.O)  A(K )=A(K ) -DOT ( A (K- I H ) ,A ( I D- I H ) ,IH ) 

K=K  +  1 

....  REDUCE  DIAGONAL  TERM 
IF(.NOT.AFAC)  GO  TO  5 
I R=JR+1 
I E  =JD  - 1 
K=J-JD 

DO  4  I  =  I R ,  IE 
ID=JDIAG(K  +  I  ) 

IF(A( ID) .EQ.O.ODO)  GO  TO  4 
D = A  ( I ) 

A( I )  =  A C I )/A( ID  ) 

A ( J  D )=A( JD )-D*A( I ) 

CONTINUE 
....  REDUCE  RHS 

I F ( BACK  )  B(J  )=8(J  ) -DO  T (A(JR+1),B{IS-1),JH-1) 

JR=JD 

IF (  .NOT  BACK)  RETURN 
....  DIVIDE  BY  DIAGONAL  PIVOTS 
DO  7  1  =  1, NEQ 
I D  = J D I  AG ( I  ) 

IF (A( ID) .ME .0.000 )  B( I )=B{ I )/A( ID) 

AENGY=AENGY+B( I )*B( I )*A( ID) 


C  ....  BACKSUBSTITUTE 

J=NEQ 

JD=JDIAG(J ) 

8  0=B( J  ) 

J=J-1 

IF(J.LE.O)  RETURN 
JR=JDIAG(J ) 

IF (JD-JR.LE . 1 )  GO  TO  10 
IS=J-JD+JR+2 
K=JR-IS+1 
00  9  I=IS,J 

9  B( I )=B( I ) -A { I+K )*D 

10  JD=JR 
GO  TO  8 
END 


SUBROUTINE  ADDSTE  ( A,B ,C , S ,P , JO IAG ,LD ,NST  ,NEL ,AFL ,BFL ,CFL  ) 
IMPLICIT  REAL*8( A-H ,0-Z ) 

....  ASSEMBLE  GLOBAL  ARRAYS 


CALLED  BY  -  FORMEL 

LOGICAL  AFL ,BFL ,CFL 
DIMENSION  A(l),  8(1),  JDIAG(l),  P( 
DO  2  J=1,NEL 
K-LD ( J  ) 

IF (K.EQ.O)  GO  TO  2 
IF(BFL)  B(K)=B(K)+P(J ) 

IF ( .NOT. AFL .AND . .NOT .CFL )  GO  TO  2 
L  =  J  D I  AG  ( K  )-K 
DO  1  1  =  1,  NEL 
M=LD(  I ) 

IF(M.GT.K.OR.M.EQ.O)  GO  TO  1 
M=L+M 

IF(AFL)  A(M)=A(I1)+S(I,J) 

IF(CFL)  C(M)=C(M)+S(J ,1 ) 

CONTINUE 

CONTINUE 


),  S(NST,1 ) ,  LD ( 1 ) ,  C(l) 


SUBROUTINE  BANADD  (GK ,MXR ,MXC ,GF ,FE ,EK ,NK ,NDFE ,LD , STI F ,FORC ) 


IMPLICIT  REAL*8( A-H ,0-Z ) 

LOGICAL  FORC.STIF 

DIMENSION  GK (MXR.MXC ) ,  GF ( 1 ) ,  FE(1),  EK(MK  ,NK) ,  LD ( 1 ) 

THIS  SUBROUTINE  ASSEMBLES  ELEMENT  ARRAYS  (EK)  AND  (FE)  TO  FORM 
GLOBAL  ARRAYS  (GK)  AND  (GF),  RESPECTIVELY.  (EK)  IS  ASSUMED  TO  BE  S 

(GK)  =  THE  UPPER  HALF  OF  A  SYMM.  GLOBAL  MATRIX  STORED  IN 

RECTANGULAR  ARRAY  WITH  DIAGONALS  ON  THE  FIRST  COLUMN. 

MXR  =  MAX.  NO.  OF  ROWS  OF  (GK)  AS  DIMENSIONED  IN  THE  MAIN. 

MXC  =  MAX.  NO.  OF  COLS  OF  (GK)  AS  DIMENSIONED  IN  THE  MAIN. 

(EK)  =  SYMMETRIC  ELEMENT  MATRIX  OF  ACTUAL  SIZE  NDFE  BY  NDFE 

SUPPLIED  BY  THE  MAIN. 


NDFE  =  DEGREE  OF  FREEDOM  FOR  THIS  ELEMENT 

=  (NO.  OF  NODES)  *  (DEGREE  OF  FREEDOM)  FOR  THIS  ELEMENT. 

(LD)  =  ELEMENT  ASSEMBLY  POINTER  ARRAY. 


BOUNDARY  CONDITIONS  WILL  BE  IMPOSED  ON  LATER  IN  THE  SUBROUTINE  BNC 

DO  2  1=1, NDFE 
I I=LD  ( I ) 

IF(FORC)  GF ( 1 1  )=GF ( 1 1 )+FE ( I ) 

IF ( .NOT. STIF  )  GO  TO  2 
DO  1  J=1 ,NDFE 
JJ=LD ( J  ) 

IF(II.GT.JJ)  GO  TO  1 
K  =J J - 1 1+1 

GK  ( 1 1  ,K  )=GK  ( 1 1  ,K  )+EK  ( I ,  J  ) 

CONTINUE 

CONTINUE 

RETURN 


M 


SUBROUTINE  BANSOL  ( A,X ,B ,N , IB .BACK ,FORW) 


IMPLICIT  P.EAL*8(A-H,0-Z) 

LOGICAL  BACK ,FORW 
DIMENSION  A( N , IB ) ,  X(l),  B ( 1 ) 

IF (BACK)  GO  TO  3 

-  THIS  SUBROUTINE  TRIANGULARIZES  A  BANDED  AND  SYMMETRIC 

COEFFICIENT  MATRIX  (A). 

-  ONLY  THE  UPPER  HALF  BAND  PORTION  OF  THE  COEFFICIENT  MATRI 

IS  STORED  AS  A  RECTANGULAR  ARRAY 
DO  2  1=2, N 
Ml=MINO( IB-1, N-I+l) 

DO  2  J=1 ,M1 
SUM=0.0 

Kl=MINO( 1-1 , IB-J ) 

DO  1  K=1,K1 

1  StJM=SUM+A( I-K ,K+1 ) *A( I-K ,J+K )/A( I-K ,1 ) 

2  A ( I , J )=A( I ,J ) — SUM 

.  THIS  SUBROUTINE  MULTIPLIES  THE  INVERSE  OF  LEFT  TRIANGULAR 

FORM  WITH  THE  RIGHT  HAND  SIDE  VECTOR,  AND  THEN  SOLVES  FOR  THE 
UNKNOWNS  BY  BACK  SUBSTITUTION  PROCESS 
IF(FORW)  RETURN 

3  NP1=N+1 

DO  5  1=2  ,N 
SUM=0.0 

Kl=MINO( IB-1 ,1-1) 

L  4  K  =  1  ,K1 

4  SUM=SUM+A( I -K ,K+1 ) /A ( I -K , 1 ) *B ( I-K ) 

5  B( I )=B( I ) -SUM 

-  BEGIN  BACK  SUBSTITUTION 

X(N)=B(N)/A(N,1) 

DO  7  K=2,N 

I=NP1-K 

J1=I+1 

J2=MINO(N , I+IB-1 ) 

SUM=0.0 
DO  6  J  = J 1  ,J2 
MM= J - J 1 +2 

6  SUM=SUM+X(J )*A( I ,MM) 

7  X ( I )=(B( I )-SUM)/A( 1,1) 

RETURN 


SUBROUTINE  BCCODE  ( 10 , IDL ,NDF ,NUMNP ,PRT  ,TRD ,TWR ,HEAD,NTAPE ) 

IMPLICIT  REAL*8(A-H,0-Z) 

LOGICAL  PRT ,TRD ,TWR 

DIMENSION  ID(NDF, 1 ) ,  IDL ( 1 ) ,  HEAD(20) 

COMMON  /LABL /  PD(9) ,A(9) ,BC(2),DUM(15 ) ,FOL (3) ,F0N(4) 

CALLS 

FOL (2 )=A(NDF+1 ) 

IF (PRT )  WRITE  (6,12)  HEAD, ( I ,BC,I=1 ,NDF) 

IF(TRD)  READ  ( NTAPE ,FOL  )  (N ,NG , ( ID ( I ,N ) , 1=1 ,NDF ) ,N=1 ,NUMNP ) 

IF(TRD)  GO  TO  6 

DO  1  1=1 ,NDF 

DO  1  J=1 ,NUMNP 

ID ( I , J  )=0 

N=0 

NG=0 

L=N 

LG=NG 

READ  (5,10)  N ,NG , ( IDL ( I ) , 1=1 ,NDF ) 

IF(N.LE.O.OR.N.GT.NUMNP)  GO  TO  6 
DO  3  1=1, NDF 
ID ( I ,N )  =  IDL ( I ) 

LG=ISIGN(LG,N-L) 

L=L+LG 

IF( (N-L )*LG.LE .0)  GO  TO  2 

DO  5  1=1, NDF 

ID( I ,L )= ID( I ,L-LG) 

GO  TO  4 
CONTINUE 

IF ( .NOT. PRT )  RETURN 
DO  9  N=1 ,NUMNP 
DO  7  1  =  1, NDF 

IF( ID ( I ,N ) . NE. 0 )  GO  TO  8 

CONTINUE 

GO  TO  9 

WRITE  (6,11)  N , ( ID( I  ,N ) , 1  =  1 ,NDF ) 

CONTINUE 

NG=0 

IF(TWR)  WRITE  (NTAPE, FOL)  (J ,NG , ( ID( I ,J ) , 1  =  1 ,NDF ) ,J  =  1 ,NUMNP ) 
RETURN 

FORMAT  (1615) 

FORMAT  (110,4113) 

FORMAT  (1H1 ,5X ,20A4,//5X , 'NODAL  B.C . 7/6X , 'NODE  1 ,9( I7,A4,A2 )/) 


A 
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SUBROUTINE  BNCOND  (S,FL , IDL  ,BVL  ,NST  ,NDFE ) 

C  . . . . . . 

IMPLICIT  REAL *8 ( A-H ,0-Z) 

DIMENSION  IDL  ( 1 ) ,  BVL(l),  S(NST,1),  FL(1) 
C 


c 

MODIFY 

THE  R.H.S.  ACCORDING  TO  ELEMENT  B.C. 

c 

BVL 

BOUNDARY  VALUE 

c 

S 

ELEMENT  STIFFNESS  MATRIX  (NST, NST) 

c 

FL 

RIGHT-HAND  SIDE 

c 

NDFE 

NUMBER  OF  EQUATIONS 

c 

NST 

MAX.  NUMBER  OF  EQUATION  PER  ELEMENT 

C 

C  CALLED  BY  -  FORMEL 

C 

DO  2  1  =  1  ,NDFE 
I BC= I DL ( I ) 

C 

C  CODE  FOR  DIRICHLET  TYPE  CONSTRAINT  BOUNDARY  CONDITIONS  ARE  GIVEN 
C  AS  '1'  OR  '-1',  AND  HAVE  BEEN  CHANGED  TO  'O'  ALREADY. 

C 

IF(IBC.NE.O)  GO  TO  2 
C 

C  MODIFY  THE  R.H.S.  FOR  NON  HOM.  DIRICHLET  COND. 

C 

DO  1  J=1 ,NDF£ 

IF(I.NE.J)  FL(J)=FL(J)-S(J,I)*BVL(I) 

1  CONTINUE 

2  CONTINUE 
RETURN 
END 


o  o  o 
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SUBROUTINE  CPTABL  (TITLE, IRITE) 

C  — - - - - 

IMPLICIT  REAL*8(A-H  ,0-Z) 

COMMON  /TIMECK/  CP  T I  ME ( 20 ) .COMPUT (20) ,NCHK 

CALLED  BY  -  ACNTRL 

NCHK=NCHK+1 

N=NCHK-1 

IF ( IRITE. EQ.O)  GO  TO  1 

WRITE  (6,4)  ( COMPUT ( I ) ,CPTIME { I ) , I =1 ,N ) 

GO  TO  3 

1  IF(NCHK.LE.20)  GO  TO  2 

WRITE  (6,4)  (COMPUT (I), CPTIME (I), 1=1, N) 

NCHK=1 

2  CPTIME (NCHK )=RCL0K1 (1.0) 

COMPUT ( NCHK )=T I TLE 

3  RETURN 
C 

4  FORMAT  (1H1.///10X, 'TIME  ELAPSED  FOR  EACH  COMMAND  7/20 (10X , 'MACRO 1 

1 ,5X ,A4 , '  ;  1  ,F10 .6 , 1  SECONDS'/)) 

END 


FUNCTION  DOT  (A,B,N) 


IMPLICIT  REAL*8 (A-H ,0-Z ) 

....  VECTOR  DOT  PRODUCT 

CALLED  BY  -  ACTCOL 

DIMENSION  A( 1 ) ,  B ( 1 ) 
DOT=O.ODO 
DO  1  1  =  1  ,N 
DOT=DOT+A( I  )*B(  I ) 

RETURN 

END 


o  o  o  o  o 


SUBROUTINE  ELOATA  ( IX , IDE  ,NEN1 ,NDF ,PRT ,ERR ,TRD ,TWR , IDIF ) 

IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON  /CDATA/  0 ,HEAD(20 ) ,NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ ,MEQ,NTAPE ,NPN 
COMMON  /LABL/  PD(9 )  ,A(9 )  ,BC(2 )  ,DI (6) ,CD(3 ) ,TE (3) ,FD(3 ) ,FOL (3 ) ,FON ( 
14) 

LOGICAL  ERR,PRT,TRD,TWR 
DIMENSION  IX ( NEN1 ,  1 ) ,  IDL(l) 

CALLED  BY  -  ACNTRL 

FOL { 2 )=A( NEN1 ) 

IF(TRD)  GO  TO  10 
L=0 

DO  7  M=1 ,NUMEL 
IF(L-N)  1,3,5 

READ  (5,15)  L,LK,LX,(IDL(K),K=1,NEN) 

L  =  ELEMENT  NO. 

LK  =  MATERIAL  SET  NO. 

LX  =  GENERATION  INCREMENT.  ALL  NODES  ARE  INCREMENTED  BY  LX 
FOR  EVERY  ELEMENT  BETWEEN  L  AND  L  ON  THE  NEXT  CARD. 
ELEMENTS  MUST  BE  IN  ORDER. 

MP1=N+1 

IF(L.EQ.O)  L=NUMEL+1 
IF(LX.EQ.O)  LX=1 
IF(L-N)  2,3,5 
WRITE  (6,18)  L,N 
ERR= .TRUE . 

GO  TO  7 
NX=LX 

DO  4  K  =  1  ,NEH1 
IX (K ,L )  =  IDL (K  ) 

I X ( NEN 1  ,L)=LK 
GO  TO  7 

IX (NEN1  ,N )  =  IX (NEN1 ,N-1 ) 

DO  6  K  =  1  ,NEN 
IX (K ,N )  =  IX (K ,N-1  )+NX 
IF(IX(K,N-l).EQ.O)  I X (K  ,N  )=0 
CONTINUE 

. COMPUTE  THE  HALF  BAND  WIDTH  IBAND . 

IDIF =0 

NEN 1 1 =NEM - 1 

00  9  I =1 ,NUMEL 

DO  9  J  =  1  ,NENI I 

IF ( I X ( J  ,  I )  .LE .0 )  GO  TO  9 

L=J  +  1 

DO  8  K  =  L  ,NEN 

IF(  IX (K , I ) ,LE .0 )  GO  TO  8 

NDIF=IA3S(IX(K,I)-IX(J,I)) 

IF ( IDIF.LT.NDIF )  I D I F=MD IF 

CONTINUE 

CONTINUE 
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IF(TWR)  WRITE  (NTAPE ,F0L )  ( J  ,  ( IX( I  ,J  } , 1  =  1 ,NEN1 } ,J=1 ,NUMEL ) , IDIF 
GO  TO  11 

10  READ  (NTAPE, FOL )  ( J  ,  ( I X ( I  ,J ) , I  =  1  ,NEN1 ) , J  =  1 ,NUMEL ) , IDIF 

11  IF(.NOT.PRT)  RETURN 
DO  12  I -1 ,NUMEL ,50 

WRITE  (6,16)  0 ,HEAD , (K ,K  =  1 ,NEN  ) 

J=MINO ( NUMEL , 1+49  ) 

DO  12  N  =  I ,  J 

12  WRITE  (6,17)  N ,  I X ( NEN 1 ,N ) , ( IX (K ,N  )  ,K  =  1  ,MEN  ) 

IF(TWR)  WRITE  (6,13)  NTAPE 

IF(TRD)  WRITE  (6,14)  NTAPE 
RETURN 
C 

13  FORMAT  ( //5X , * . SAVED  ON  FILE  NO  : ' , 14) 

14  FORMAT  (//5X,‘ . RTEAD  FROM  FILE  NO  : ' , 14 ) 

15  FORMAT  (1615) 

16  FORMAT  ( 1H1 ,A1 ,5X ,20A4//5X .8HELEMENTS//3X .7HELEMENT ,2X .8HMATERI AL , 
114(13, 5H  NODE  )/( 20X , 14( 13 ,5H  NODE))) 

17  FORMAT  (2110, 14I8/(20X, 1418)) 

18  FORMAT  ( 5X ,20H**ERR0R  03**  ELEMENT , 15 ,22H  APPEARS  AFTER  ELEMENT, 15 
1) 

END 


SUBROUTINE  ELMLIB  ( IXL ,XL ,D ,FL ,S ,UL ,TL ,NDM,NDF ,NST ,NEL ,ISW,N  ,IEL ,N 


1G) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  IXL ( 1 ) ,  XL ( 1 ) ,  D(l),  FL(1),  S(NST,1),  UL(1),  TL { 1 ) 

CALLED  BY  -  FORMEL 
CALLS  -  ZERO 


IF( ISW.EQ .2 )  CALL  ZERO  (S,NST*NST) 

IF( ISW.NE . 5 )  CALL  ZERO  (FL.NST) 

IF ( IEL -EQ.2  )  GO  TO  2 
IF( IEL.EQ.3)  GO  TO  3 

CALL  ELMT01  ( IXL ,XL ,D ,FL ,S ,UL ,TL ,NDM,NDF ,NST ,NEL , ISW, IEL ,N ,NG ) 
RETURN 

CALL  ELMT02  ( IXL ,XL ,D ,FL ,S ,UL ,TL .NDM.NDF ,NST ,NEL , ISW, IEL ,N ,NG ) 
RETURN 

CALL  ELMT03  ( IXL , XL ,D,FL,S, UL.TL, NDM.NDF, NST.NEL, ISW, IEL, N,NG) 
RETURN 
END 


SUBROUTINE  FORMEL  (X ,XL ,D ,DL ,S ,UL ,FL ,TL ,U ,F ,T , I X , IE , ID ,LD ,  I  XL ,  JDIA 


1G,NDF ,NM,NEN1 ,NST , ISW,A,C,B ,AFL ,BFL ,CFL ,PRT ,TLD,NRD  , 
2NOW,IBAND,BND) 

IMPLICIT  REAL*8{A-H ,0-Z ) 

LOGICAL  AFL ,BFL ,CFL ,DFL ,EFL ,PRT ,TLD,BND 

COMMON  /CDATA/  0 , HEAD (20 ) .NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ ,MEQ ,NTAPE ,NPN 
COMMON  /STRES/  EPS (6 )  ,SIG{6) 

DIMENSION  X (NM, 1 ) ,  D(NRD .NUMMAT  ,1 ) ,  F(NDF,1),  XL(NM,1),  FL(NDF.l), 
1  S ( NST , 1 ) ,  T ( 1 ) ,  UL(NDF.l),  TL ( 1 ) ,  IX(NENl.l),  LD(NDF.l),  ID(NDF,1 
2),  JDIAG(l),  B ( 1 ) ,  A( 1 ) ,  C(l),  U(l),  IXL(l),  IE(4,1),  DL(NRD) 

CAI LED  BY  -  ACNTRL 

CALLS  -  ELML I B .BNCOND , ADDSTF , BANNAD 

NDM=NM 

MA1=0 

IF( ISW.EQ.5 )  REWIND  14 
REWIND  13 

....  LOOP  ON  ELEMENTS 
DO  11  N=1 .NUMEL 
MA= I X ( NEN 1 , N ) 

IF(MA.EQ.MAl)  GO  TO  2 

I  EL  =  IE ( 1 ,MA) 

NT=IE (2 ,MA) 

NGP=IE(3,MA) 

IF( ISW.EQ.5)  NGP=IE(4,MA) 

DO  1  1  =  1, NRD 
DL ( I )=D ( I ,MA , 1 ) 

MA1=MA 
EFL=. FALSE. 

....  SET  UP  LOCAL  ARRAYS 
DO  8  1=1, NEN 

I I  =  I X  ( I  ,N) 

I  XL ( I )  =  1 1 
TL( I )=0.0D0 
IF(II.NE.O)  GO  TO  5 
DO  3  J=1,NDM 

XL ( J , I )=0.0D0 
DO  4  J  =  1  ,NDF 
UL(J,I)=O.ODO 
LD( J , I )=0 
GO  TO  8 
NEL  =  I 

I  ID=I I*NDF-NDF 

IF{ TLD.OR. ISW.EQ .2 )  TL ( I  )=T( 1 1 ) 

00  6  J  =  1  ,NDM 
XL ( J  ,  I )=X ( J ,11) 

DO  7  J=1,NDF 
K  =  IABS ( ID ( J , 1 1  ) ) 

UL( J  , I  )=F( J ,11) 

IF (K.GT.O)  UL (J  , I  )=U(K  ) 
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IF (K.EQ.O.AND.F (J , 1 1 J.NE.O)  EFL= .TRUE . 

LD(J,I  )=K 

IF(BND)  UL ( J , I )=U ( 1 1 0+J  ) 

7  IF(BND)  LD ( J , I )  =  1 1 D+J 

8  CONTINUE 
C 

CALL  ELMLIB  ( I  XL  ,XL  ,DL  ,FL  ,S  ,IJL  ,TL  ,NDM,NDF  ,NST  ,NEL  ,  ISW,N ,  IEL ,NGP ) 

C 

NDFE=NDF*MEL 

IF ( . NOT.PRT )  GO  TO  9 

WRITE  (6,14)  N,((LD(I,J),I=1 ,NOF ) ,J=1 ,NEL) 

IF ( ISW.EQ.2 )  WRITE  (6,15)  ( (S( I  ,J ) ,J=1 ,NST ) , 1  =  1 ,NST ) 

I F ( ISW.EQ.5)  WRITE  (6,16)  ( ( UL ( I ,J ) , 1  =  1 ,NDF )  ,J  =  1  ,NEL ) 

IF ( ISW.EQ.3.0R. ISW.EQ.4)  WRITE  (6,17)  ( (FL ( I  ,J ) , 1=1 ,NDF ) ,  J  =  1 , NEL  ) 
IF( ISW.EQ. 3  )  WRITE  (6,18)  (TL ( I ) ,  1  =  1  ,NEL ) 

9  IF(BND)  GO  TO  10 

IF (AFL . AHD.BFL .AND .EFL  )  CALL  BNCOND  (S.FL ,LD ,UL ,NST,NDFE ) 
IF(AFL.0R.3FL .OR.CFL  )  CALL  ADDSTF  ( A,B ,C , S ,FL , JDI AG ,LD ,NST ,NDFE , AF 
1L ,BFL ,CFL  ) 

GO  TO  11 

10  IF (AFL -OR.BFL )  CALL  BANADD  ( A ,MEQ ,  IBAND ,B ,FL ,S ,NST ,NDFE  ,LD , AFL ,BFL 
1) 

11  CONTINUE 

IFdSW.NE. 5. OR.  .NOT.PRT)  GO  TO  13 
IF ( IEL .LE.2 )  GO  TO  13 
REWIND  14 

DO  12  1 1 =1 ,NUMEL ,50 
WRITE  (6,26)  HEAD 
JJ=MINO(NUMEL , 11+49 ) 

DO  12  J  =  1 1  ,JJ 
MA=IX (NENl  ,J  ) 

NG=IE (4,MA) 

I EL  =  I E ( 1 ,MA) 

DO  12  K  =  1  ,NG 

READ  (14)  XX,YY,ZZ,EPS,P1,P2,P3,A1 

12  CONTINUE 

13  CONTINUE 

I F (  .NOT.PRT. OR. I SW.NE.2)  RETURN 
JKL  =JD I  AG ( NEQ ) 

IF(BND)  jkl=meq*iband 

WRITE  (6,19)  ISW,NEQ,MEQ,JKL 

IF(.NOT.BND)  WRITE  (6,20)  (JDI AG ( I  )  ,1  =  1  ,NEQ) 

WRITE  (6,21)  (A( I ) , I  =  1 ,JKL  ) 

IF(CFL)  WRITE  (6,22)  (C ( I  ) .  1  =  1 ,  JKL  ) 

WRITE  (6,23)  ( B ( I ) , I =1 ,NEQ ) 

WRITE  (6,24)  (U ( I ) , I =1  ,NEQ  ) 

WRITE  (6,25)  ( (F( I ,J ) , 1=1 ,NDF ) ,J=1 ,NUMNP ) 

RETURN 

C 

14  FORMAT  (//2X,  'STIFF.  N  =  ',I5,5X,‘LD  ‘,2414//) 

15  FORMAT  ( 10X , 8E 12.4) 

16  FORMAT  (5X,’IJL',3X,8E12.4) 


17  FORMAT  ( 5X , ' FL 1 , 3X ,8E 12 .4  ) 

18  FORMAT  (5X,'TL‘ ,3X,8E12.4) 

19  FORMAT  (//SX.’ISW,  NEQ ,  MEQ ,  JKL  =  NCOEF , ' ,4 1 10/ ) 

20  FORMAT  ( /5X , 1  JO IAG  1  ,2016  ) 

21  FORMAT  (5X , 1 A  1 ,2X , 10E12 .4 ) 

22  FORMAT  (5X,'C  ' ,2X .10E12.4) 

23  FORMAT  (5X,'B  ' ,2X , 10E12.4) 

24  FORMAT  (5X,'U  ' .2X.10E12.4) 

25  FORMAT  (5X,'F  ' ,2X , 10E12 . 4 ) 

26  FORMAT  ( 1H1 , 5X ,20A4 , / , 5X  , ' _ ELEMENT  STRAINS  _  CCW  ROTATION 

r  IS  POSITIVE  _ 1 ,/,5X, ' _  1-DIRECTION  IS  RADIAL  DIRECTION,' 

2  2-DIRECTION  IS  AXIAL  DIRECTION  FOR  AXISYMMETRIC  PROBLEMS  _ 1 , 

3/,2X , ' EL  1-COORD.  2-COORD.  1-STRAIN  2-STRAIN  3-STR ' , ' AI 
4N  12-STRAIN ' ,T78,'P1‘ ,T89,'P2' ,T100,'P3' .Till,' Al * ,/) 


SUBROUTINE  GENVEC  (JKL  ,X ,CDO  ,PRT .ERR ,TRD ,TWR ) 


IMPLICIT  REAL*8(A-H ,0-Z) 

....  GENERATE  REAL  OATA  ARRAYS  BY  LINEAR  INTERPOLATION 

CALLED  BY  -  ACNTRL 
CALLS  -  PCOMP ,  ZERO 
LOGICAL  PRT, ERR, PCOMP, TRD.TWR 

COMMON  /CDATA/  0 .HEAD (20 ) ,NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ ,MEQ ,NTAPE ,NPN 
COMMON  /L ABL /  PD (9 )  ,A(9 )  ,BC (2 ) ,DI (6 ) ,CD(3 ) ,TE (3 ) ,FD (3 ) ,FOL ( 3 ) ,FON ( 
4) 

DIMENSION  X ( JKL  ,  1 ) ,  XL (6 ) ,  CDD (3 ) ,  FACT(6) 

DATA  BL/4UBLAN/,GE0M/4H  COO/ 

FON ( 2 ) = A ( JKL ) 

IF(TRD)  GO  TO  12 
CALL  ZERO  ( X  , JKL*NIJMNP  ) 

IF(. NOT. PCOMP (CDD (1),GE0M))  GO  TO  2 

DO  1  N=1 ,NUMNP 

X ( 1 , N )  =  BL 

N=0 

NG=0 

L=N 

LG=NG 

IF ( LG  .EQ.  0)  READ ( 5 ,1000 )  N ,NG ,  (XL ( I ) , I =1 ,JKL ) ,  (FACT (J ) , J  =  1 ,JKL ) 
IF(LG.EQ.O)  READ  (5,17)  N ,NG ,  (XL ( I ) ,  1=1 ,  JKL  ) 

IF(LG.NE.O)  READ  (5,17)  N ,NG , (XL ( I ) ,  1  =  1 ,JKL  ) 
IF(M.LE.O.OR.N.GT.NUMNP)  GO  TO  11 
DO  4  1  =  1, JKL 
X ( I ,M  )=XL  ( I ) 

IF (LG )  5,3,5 
LG  =  I S I GN (LG ,N-L ) 

L I  =  ( I ABS (N-L+LG )-l ) /I ABS (LG ) 

DO  7  1  =  1, JKL 
F=1 . +FACT ( I  ) 

F0=1 .000 
FF=FO 

DO  6  J=2  ,  L I 

FO=FO*F 

FF=FF+FO 

XL(I)=(X(I,N)-X(I,L) )/FF 

FF =1.000 

L=L+LG 

IF ( (N-L )*LG.LE.O)  GO  TO  3 

IF (L .LE .O.OR.L.GT. NUMNP )  GO  TO  10 

DO  9  1=1, JKL 

X ( I ,L )  =  X ( I ,L-LG  )+XL ( I  ) 

XL(I)=XL(I)*(1. +FACT ( I  ) ) 

GO  TO  8 

WRITE  (6,21)  L , (CDD ( I ) , 1  =  1 ,3 ) 

ERR=.TRUL. 


GO  TO  3 

IF(TWR)  WRITE  ( NTAPE  ,FON  )  ( J , (X ( I ,J ) ,1  =  1 , JKL ) , J  =  1 , NUMNP ) 

GO  TO  13 

READ  (NTAPE, FON)  ( J , (X ( I ,J ) , 1=1 ,JKL ) ,J=1 .NUMNP ) 

IF ( .NOT.PRT )  RETURN 
DO  14  1=1, NUMNP, 50 

WRITE  (6,18)  HEAD, (CDD(L ) ,L=1 ,3 ) , (L ,CDD(1 ) ,CDD (2 ) ,L=1  ,JKL ) 

N=MI NO (NUMNP ,1+49 ) 

DO  14  J= I ,N 

I F ( PCOMP ( X ( 1 , J ) ,BL ) )  WRITE  (6,19)  N 

IF( .NOT.PCOMP (X ( 1 ,J ) ,BL ) )  WRITE  (6,20)  J , (X (L  , J  )  ,L  =  1  ,JKL ) 

IF(TWR)  WRITE  (6,15)  NTAPE 
IF(TRD)  WRITE  (6,16)  NTAPE 
RETURN 

FORMAT  (//5X, ' . SAVED  ON  FILE  NO  :',I4) 

FORMAT  ( //5X , ' . READ  FROM  FILE  NO  :',I4) 

FORMAT  (2I5.7F10.0) 

FORMAT  ( 1H1 ,5X .20A4//5X .5HN0DAL .3A4//6X .4HN0DE ,9( 17 ,A4,A2 ) ) 

FORMAT  (I10.32H  HAS  NOT  BEEN  INPUT  OR  GRNERATED ) 

FORMAT  ( 110, 9F 13.4) 

FORMAT  ( 5X ,43H**FATAL  ERROR  02**  ATTEMPT  TO  GENERATE  N0DE,I5,3H  IN 


o  o  o  o  c~> 
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SUBROUTINE  MATLIN  (D ,NADD ,IE .NUMMAT ,NRD,DL ) 


....  CALLEO  BY  :  ACNTRL 
....  CALLS  :  MATRLD,  MELAW 
READ  ALL  MATERIAL  DATA  AND  STORE  IN  (D). 

FORM  STRESS  -  STRAIN  LAW  IF  ISOTHERMAL  MATERIAL  (NTEMP  =1). 
IMPLICIT  REAL*8(A-H ,0-Z ) 

REAL *8  K1,K2 
LOGICAL  STRES 

COMMON  /ML I ST/  IEL , NTEMP .NGAUSS .NSTRES ,AHG ,T , El ,E2 ,XNU1 ,XNU2 ,G2 ,AL 
1PA1 ,ALPA2 , SC, ST, PHI ,RHO ,K1 ,K2 ,C1 ,C2 ,A,C,A1 ,A2,A3 ,A4 
COMMON  /ELABL/  I  TYPE (6, 10) 

DIMENSION  D(NRD , NUMMAT ,1 ) ,  IE(4,1),  MHEAD(20) ,  DL(NRD) 

NAMELIST  /MATRL/  T,E1 ,E2 ,XNU1 ,XNU2,G2,RH0,ALPA1 ,ALPA2 ,SC , ST, PHI ,K1 

NAMELIST  /MATID/' MATNO ^ELTYPE .NTEMP .NGAUSS .NSTRES ,ANG 
MTEMP=1 
WRITE  (6,7) 

C 

00  3  MA=1, NUMMAT 

C 

READ  (5,6)  MHEAD 
READ  (5.MATID) 

IEL=ELTYPE 
STRES=. FALSE. 

STRES= .TRUE . 

MTEMP=MAXO(MTEMP , NTEMP ) 

IF (MATNO.LT. l.OR.MATNO.GT. NUMMAT)  GO  TO  4 

IE ( 1 ,MATNO)=IEL 

IE(2,MATN0)=NTEMP 

IE ( 3 ,MATNO ) =NGAUSS 

IE (4 ,MATNO)=HSTRES 

C  IF ( STRES )  WRITE (6 ,610 )  MATNO ,MHEAD , IEL , ( ITYPE ( I ,IEL ) , 1=1 ,6 ) 

IF(STRES)  WRITE  (6,8)  MATNO, MHEAD, IEL 
IF(. NOT. STRES)  WRITE  (6,8)  MATNO .MHEAD, IEL 
WRITE  (6,9)  NGAUSS, NSTRES, NTEMP, ANG 
C 

DO  2  NT=1, NTEMP 

G2=0.0D0 

READ  (5, MATRL) 

IF(G2.LE.0.0D0)  G2=0.5D0*El/( 1 .0D0+XNU1 ) 

IF(STRES)  GO  TO  1 
C 

C  ....  MATERIAL  PROPERTIES  FOR  A  TRANSPORT  PROBLEM 
CALL  MATRLD  (14,1, 6, MATNO, NT ,D, NUMMAT, NRD) 

C 

D(7,  MATNO,  NT)=RHO 

WRITE  (6,10)  (D(I, MATNO, NT),  1=1, 7) 

C 

1  CALL  MATRLD  (1,11,16 , MATNO , NT ,D,NURMAT ,NRD ) 

WRITE  (6,11)  NT, (D(I, MATNO, NT), 1=12, 26) 

C 
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2  CONTINUE 
C 

IF( .NOT.STRES.OR.NTEMP .GE.2 )  GO  TO  3 

C  _  ISOTHERMAL  MATERIALS  -  FORM  ELASTIC  LAW  HERE 

CALL  MELAW  ( IEL ,D( 1 , MATNO ,1 )  ,DL  ) 

3  CONTINUE 
C 

NADD=NUMMAT*MTEMP*NRD 

C 

GO  TO  5 

4  WRITE  (6,12)  MATNO 
STOP 

5  RETURN 
C 

6  FORMAT  (20A4) 

7  FORMAT  (///5X,'M  ATERIAL  DATA'/) 

8  FORMAT  (//5X , 'MATERIAL  SET  NO.  =  * , 14 ,  *  :  ' .20A4//5X , 'ELEMEN 

IT  TYPE  NO.  =',I4,'  :  ',6A4) 

9  FORMAT  (5X,'NO.  OF  GAUSS  PTS.  = ' ,I4/5X , 'NO.  OF  STRESS  PTS. 

1  = ' , I4/5X , ' NO .  OF  TEMP  -  MATRL  PTS. =' .I4/5X , 'ANGLE  FROM  XI  TO  E 

21  = ' ,F10.3 , '  OEGREE  (  C.  C.  W.  IS  POSITIVE).'/) 

10  FORMAT  (8X, 'CONDUCTIVITY  K1  =  ' .E14.5/8X , 'CONDUCTIVITY  K2 

1  =' .E14.5/8X, 'CONVECTION  Cl  =' .E14.5/8X, 'CONVECTION  C 

22  =' .E14.5/8X, 'SOURCE/SINK  A  =' .E14.5/8X , 'CAPACITY 

3  C  =',E14.5/8X, 'MASS  DENSITY  RH0=',E14.5) 

11  FORMAT  (5X , 'TEMP.  PT.  NO.  =',I4,'  :  TEMPERATURE  =' .E14.5//8X, 'MOD. 

1  OF  ELASTICITY  El  =' ,E14.5/8X, 'MOD.  OF  ELASTICITY 

2  E2  =',E14.5/8X, 'POISSONS  RATIO  XNU1  =  ' ,E14.5/8X, 'P 

30ISS0NS  RATIO  XNU2  ='  .E14.5/8X, ' INDEP.  SHEAR  MOD. 

4  G2  =',E14.5/8X, 'FLUID  COMPRESSIBILITY  ALPA1  =',E14.5/8X 

5, 'SOLID  COMPRESSIBILITY  ALPA2  =' .E14.5/8X ,' STIFFNESS  DAMPING  C 
60EF .  SC  =' .E14.5/8X, 'MASS  DAMPING  COEF.  ST  =',E14.5 

7/8X , 'PROSOITY  PHI  ='  .E14.5/8X, 'MASS  DENSITY 

8  RHO  =  ' ,E14.5/8X, 'PERMEABILITY  K1  =',E1 

94. 5/8X, 'PERMEABILITY  K2  =' .E14.5/8X ,' FLUID  DENSITY 

£  Cl  = ' ,E14.5/) 

12  FORMAT  (//// 5X,' .  FATAL  ERROR  IN  THE  NAMELIST  INPUT  *MATID*  .. 

1...710X,  'MATERIAL  SET  NO.  (MATNO)  = ' ,  15 , '  IS  NOT  POSSIBLE. 

2  EXECUTION  TERMINATED. ' ) 

END 


ooooooooooooooooooo 
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SUBROUTINE  MATRLD  (NCI ,ND1 ,N ,MA,NT,D ,NUMMAT,NRD) 

C  . . - . - . — . 

IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON  /MUST/  NDUM(4 )  ,DUMV  (23 ) 

DIMENSION  D(NRD,NUMMAT,1) 

....  CALLED  BY  :  MATLIN 

PUT  NAMELIST  INPUT  VARIABLES  (  STORED  IN  DUMY  )  INTO  THE  ARRAY  D 
D ( 1 1 ,MA ,NT )  =  ANG 
D ( 12 ,MA ,NT  )  =  T 
D ( 1 3 , MA , NT )  -  El 
D ( 14 ,MA ,NT )  =  E2 
0(15,  MA, NT)  =  XNU1 
D ( 16 ,MA ,NT  )  =  XNU2 
D ( 17 ,MA ,NT )  =  G2 

D(18,MA,NT)  =  ALPA1  :  FLUID  COMPRESSIBILITY 

D(19,MA,NT)  =  ALPA2  :  SOLID  COMPRESSIBILITY 

D ( 20 ,MA ,NT )  =  SC  :  STIFFNESS  DAMPING  COEFFICIENT 

D(21,MA,NT)  =  ST  :  MASS  DAMPING  COEFFICIENT 

D (22 ,MA ,NT )  =  PHI  :  POROSITY 

D (23 ,MA,NT )  =  RHO 

D(24,MA,NT )  =  K1  :  PERMEABILITY  IN  1 
D (25  ,MA ,NT )  -  K2  :  PERMEABILITY  IN  2 
D ( 26 ,MA,NT )  =  Cl  :  FLUID  DENSITY 

DO  1  1  =  1, N 
II=I+NC1-1 
ID=I+ND1-1 

1  D ( ID ,MA ,NT )=DUMY (II) 

RETURN 

END 


SUBROUTINE  MELAW  (IEL.D.DL) 


IMPLICIT  REAL*8(A-H,0-Z) 

....  CALLED  BY  :  MATLIN 
....  CALLS  :  METRAN 

COMPUTE  STRESS-STRAIN  RELATIONS  USING  CONSTANTS  IN  D ( 11-23 ) 

DIMENSION  D(l),  DL(1) 

E1=D(13) 

E2=D ( 14 ) 

XNU1=D(15 ) 

XNU2=D(16) 

XN=E1/E2 

1=2 

D(8)=D( 18)*(1 ,ODO+XNUl*( 1-1 ) ) 

0(9 )=D(19 )*(1 .ODO+XNU2*( 1-1 ) ) 

DD=E2/( 1 .0D0+ ( 1-1 )*XNUl )/( l.ODO- ( 1-1 )*XNU1-I*XN*XNU2*XNU2 ) 

D ( 1 )=XN*(1 .ODO- ( 1-1 )*XN*XNU2*XNU2 )*DD 
0(2 )=XM*XNU2*( 1 ,0D0+( 1-1 )*XNU1 )*DD 
0(3 )= (1 .ODO- ( 1-1 )*XNU1*XNU1 )*DD 
D ( 4 )  =0 (17) 

D(5)=0.0D0 
0 (6  )=0.0D0 
D(10  )=0.0D0 

WRITE (6, 600)  (D( I ) ,  1  =  1 ,25  ) 

RETURN 


LOGICAL  FUNCTIONPCOMP ( A,B ) 
IMPLICIT  REAL*8(A-H  ,0-Z ) 


CALLED  BY  -  GENVEC , ACNTRL ,PRTDIS 
CALLS 

PCOMP=. FALSE. 

....  IT  MAY  BE  NECESSARY  TO  REPLACE  THE  FOLLOWING  ALPHANUMERIC 

_  COMPARISON  STATEMENT  IF  COMPUTER  PRODUCES  AN  OVERFLOW 

IF (A.EO.B )  PCOMP  =  .TRUE  . 

RETURN 

END 


o  o  o  o  o  o 
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SUBROUTINE  PLOAD  (ID,F,B,R,T,NN,P,IW,BND) 


IMPLICIT  REAL*8(A-H ,0-Z ) 
LOGICAL  BND 


....  FORM  LOAD  VECTOR  IN  COMPACT  FORM 


CALLED  BY  -  ACNTRL 
CALLS  -  ZERO 


P 


DIMENSION  I D ( 1 ) ,  F(l),  B(l),  T(l),  R(l) 

IF (BND )  GO  TO  4 
IF(IW.NE.O)  GO  TO  2 

C  _ CLEAR  THE  R.H.S.(B)  AND  GET  NODAL  BOUNDARY  LOADING. 

CALL  ZERO  (B ,NN  ) 

DO  1  N=1 ,NN 
J  =  I D  ( N ) 

1  IF(J.GT.O)  B(J )=F(N)*P+R(J ) 

RETURN 

C  _ STORE  PRESCRIBED  TEMPERATURE  AND  SOLUTION  IN  (T). 

2  DO  3  N=1 ,NN 
J = I D  ( N ) 

T(N)=F(N) 

3  IF(J.GT.O)  T(N)=B{J ) 

RETURN 

4  I F ( IW.NE.O)  GO  TO  6 
CALL  ZERO  (B.NN) 

DO  5  1=1, NN 

5  IF(ID(I).EQ.O)  B( I )=F ( I )*P+R( I ) 

RETURN 

6  00  7  1=1, NN 

7  T ( I )=B( I ) 

RETURN 

END 
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SUBROUTINE  PROFIL  («JD  IAG ,  ID  ,IX  ,NDF  ,NENI  ,NAD ,TRD,TWR ) 
IMPLICIT  REAL*8(A-H ,0-Z) 


....  COMPUTE  PROFILE  OF  GLOBAL  ARRAYS 

COMMON  /CDATA/  0 ,HEAD(20 ) , NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ ,MEQ,NTAPE ,NPN 
LOGICAL  TRD,TWR 

DIMENSION  JDIAG(l),  ID(NDF  ,1 ) ,  IX(NEN1,1) 

....  SET  UP  THE  EQUATION  NUMBERS 

CALLED  BY  -  ACNTRL 

MEQ=0 

DO  4  N=l, NUMNP 
DO  3  1=1 ,NDF 
J  =  ID( I  ,N ) 

IF ( J  )  2,1,2 
NEQ=NEQ+1 
ID ( I ,N  )=NEQ 
JDIAG(NEQ)=0 
GO  TO  3 
ID ( I  ,N  )=0 
CONTINUE 
CONTINUE 

....  COMPUTE  COLUMN  HEIGHTS 
DO  9  N=1 ,NUMEL 
DO  8  1  =  1, NEN 
1 1  =  I X ( I ,N) 

IF(II.EQ.O)  GO  TO  8 
DO  7  K=1,NDF 
KK=ID(K,II ) 

IF(KK.EO.O)  GO  TO  7 
DO  6  J  =  I  ,MEN 
JJ  =  IX ( J  ,N  ) 

IF( JJ.EQ.O)  GO  TO  6 
DO  5  L=1 ,NDF 
LL  =  ID(L  ,JJ ) 

IF(LL.EQ.O)  GO  TO  5 
M=MAXO(KK  ,LL  ) 

JDIAG(M)=MAXO(JDIAG(M) , IABS(KK-LL ) ) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

....  COMPUTE  DIAGONAL  POINTERS  FOR  PROFILE 
NAD=1 

JDIAGC1 )=1 
IF(NEQ.EQ.l)  RETURN 
DO  10  N=2  ,NEQ 

JUI AG(N )=JDIAG(N )+JD IAG ( M - 1 ) + 1 
NAD=JDIAG(NEQ ) 


o  o  o  o  o  o 
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SUBROUTINE  PRTOIS  ( ID ,X ,U ,F,NM,NF ,TRD ,TWR ,PRT ,NSTEP ,  INCR .TIME ,BND ) 

C  - - - - 

IMPLICIT  REAL*8 ( A-H ,0-Z ) 

....  OUTPUT  NODAL  VALUES 

CALLED  BY  -  ACNTRL , CONSOL 
CALLS  -  PCOMP 

COMMON  /CDATA/  0 , HEAD (20) ,NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ  ,MEO  ,NTAPE  ,NPN 
COMMON  /LABL/  PD(9 ) ,A(9 )  ,BC(2  )  ,DI (6 ) ,CD(3 ) ,TE ( 3 )  ,FD(3 ) ,FOL (3 ) ,FON( 
14) 

COMMON  /PRLOD/  PROP 
LOGICAL  PCOMP ,TRD, TWR, PRT.BND 

DIMENSION  X(NM,1 ) ,  F(NF,1),  ID(NF.l),  U(l),  UL(6),  BV ( 2 ) ,  I BC ( 6  ) 
DATA  BL/4HBLAN/.BV/4H  B.V.4HALUE/ 

C 

NTAPE=1 1 
NDM=NM 
PD (2 ) = A ( NF ) 

PD  (4  )=A(!JM) 

PD (f )=A(MF ) 

PD(£)=A(NF ) 

C 

IF(TWR)  WRITE  (NTAPE ,4 )  NF 
DO  3  1 1  =  1 ,NUMNP ,50 

IF(PRT)  WRITE  (6,5)  HEAD ,NSTEP , INCR .TIME ,BC ( 1 ) ,BC (2 ) , ( I ,CD ( 1 )  .CD (2 
1 ) ,1=1 ,MDM) , ( I ,BV (1 ) ,BV(2 ) , 1=1 ,NF) , ( I ,DI (1 ) ,DI (2 ) , 1=1 ,NF) 

C  IF ( 1 1  .EQ.  1  .AND.  TWR)  WRITE(NTAPE ,610)  HEAD ,NSTEP , INCR .TIME 

IF ( 1 1 .EQ . 1 . AND .TWR )  WRITE  (NTAPE, 5)  HEAD, NSTEP, INCR, TIME, BC(1 ) ,BC( 
12),(I,CD(1),CD(2),I=1,NDM),(I,BV(1),BV(2),I=1,NF),(I,DI(1),DI(2),I 
2=1, NF) 

JJ=MINO(NUMNP, 11+49) 

DO  3  N  =  1 1  ,JJ 

IF (PCOMP (X ( 1 ,N ) ,BL ) )  GO  TO  2 
DO  1  1=1, NF 
UL ( I )=F( I ,N )*PROP 
K  =  I ABS ( ID( I ,N) ) 

I BC ( I )=0 

IF(K.EQ.O)  IBC ( I  )  =  1 
IF(.NOT.BND)  GO  TO  1 
IBC ( I  )=K 
K  =  ( N  —  1 )*NF+I 

1  IF(K.GT.O)  UL ( I  )=U(K) 

IF(PRT)  WRITE  (6, PD)  N , ( IBC ( I ) , 1=1 ,NF ) , (X ( I ,N ) , I =1 ,NDM) , (F( I ,N ) , I  = 
11  ,NF ) ,  (UL ( I ) , 1  =  1  ,NF  ) 

IF(TWR)  WRITE  (NTAPE, PD)  N , { IBC ( I ) , I =1 ,NF ) , (X ( I ,M ) , I =1 ,NDM) , (F( I ,N 
1 ) , I =1 ,NF ) ,  (UL ( I ) ,  1  =  1 ,NF ) 

2  CONTINUE 

3  CONTINUE 

C  .  '  TIME  =' .E12.4//1X, 'NODE  1 , IX, A4.A2, 12, A4.A2, 9(14, A4,A2) ) 

C  610  FORMAT ( 20A4 ,3X , 1  STEP  NO.  =',I5,‘  INCR.  NO.  =  ' ,  15 , '  TIME  =',E1 


« 


■*. 


1 


OOO  O  OC“>OC->OOOC-)OC~> 
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SUBROUTINE  PVALUE  ( S ,  IS  , N  ,X , Y  ,Z  ,LN0 ,  IEL  ,TE  ,D  ) 

C  - 

IMPLICIT  REAL*8(A-H ,0-Z ) 

C 

COMMON  /CDATA/  0, HEAD (20 ) .NUMNP ,NUMEL .NUMMAT ,NEN ,NEQ ,MEQ ,NTAPE  ,NPN 
COMMON  /ANGLE/  ANG (4 ) ,DC ( 3 ) ,ST (2 ) , IEPS.LINE ,TWR 
DIMENSION  S (6 ) ,  0(1),  IFAIL (4 ) 

DATA  RT2/1. 41421356237309/, PI23/2. 09439510239321/ 

DATA  IFAIL/1H  , 1HT , 1HC , 1HS/ 

LOGICAL  TWR 

C  ....  COMPUTE  AND  PRINT  PRINCIPAL  VALUES  OF  STRESS  (IS=1)  OR 
C  STRAIN  ( IS=2 ) 

S  =  STRAIN  OR  STRESS  TENSOR  STORED  IN  VECTOR  FORM 
IS  =  1  IF  S  IS  STRESS 
=  2  IF  S  IS  STRAIN 

N  =  ORDER  OF  THE  TENSOR  S  (  2  OR  3  ) 

....  COMPONENTS  OF  S  MUST  3E  IN  THE  FOLLOWING  ORDERS 

S 1 1 , S 1 2 ,S22 ,S13 ,S23  ,S33 

IF( IS.EQ.2)  GO  TO  2 

IF(LINE.NE.O)  GO  TO  1 
WRITE  (6,7)  HEAD 

I F ( N . EQ . 2 )  WRITE  (6,6)  ST (IS) ,ST (IS  )  ,ST ( IS  ) 

IF(N.EQ.3)  WRITE  (6,8)  ST( IS  )  ,ST( IS  )  ,ST( IS )  ,ST ( IS ) 

1  L INE=L INE+1 

2  IF(N.EQ.3)  GO  TO  3 

_  2ND  ORDER  TENSOR  _ 

S12=S(2)/IS 
SS=(S( 1 ) -S ( 3 ) )/2.0D0 
TT=(S( 1 )+S(3 ) ) /2.0D0 
RR=DSQRT(SS*SS+S12*S12 ) 

P1=TT+RR 
P2=TT-RR 
Al=45 .ODO 

IF (S12.EQ.O.ODO.OR.SS.NE.O.ODO)  A1=DATAN2(S12 ,SS)*ANG(1 )/2.0D0 

IF (LINE.EQ.50)  LINE=0 

IF(IS.NE.l)  GO  TO  5 

P3=0.0D0 

GO  TO  5 


_  3RD  ORDER  TENSOR  _ 

3  RR=O.ODO 

U=(S(1 )+S(3)+S(6) )/3.0D0 

V=S(1)*(S(3)+S(6))+S(3)*S(6)-S(2)*S(2)-S(4)*S(4)-S(5)*S(5) 

W=S(1)*S(3)*S(6)+2.0D0*S(2)*S(4)*S(5)-S(1)*S(5)*S(5)-S(3)*S(4)*S(4 

1)-S(6)*S(2)*S(2) 

T  T= 3 . 0D0*U*U-V 
SS=0 .ODO 


IF (TT.EQ. 0.000)  GO  TO  4 
UU=DSQRT(2.0D0*TT/3.0D0) 

AA=  ( W+  ( TT-U*ll  )*U )  *RT2/UU**3 
RR=DSQRT ( DABS ( 1 . 0D0-AA*AA ) ) 

RR=DATAN2(RR,AA)/3.0D0 
P1=U+U'J*RT2*DC0S(RR ) 

P2=U+UU*RT2*DC0S(RR-PI23 ) 

P3=U+UU*RT2*DC0S(RR+PI23 ) 

IF(LINE.GE.50)  LINE=0 

IF(IS.EQ.l)  WRITE  (6,9)  LNO , I  FA  I L ( IM) ,X ,Z ,S(1 ) ,S(6 ) ,S (3 ) ,S(4 ) ,P1 ,P 
12  ,P3 ,A1 ,D(12 ) ,0(13) 

IF(TWR.AND.IS.EQ.l)  WRITE  (MTAPE.10)  LNO.IFAIL(IM) ,X,Y,Z,P1 ,P2,P3, 
1A1 ,A2 ,A3 ,TE  ,D ( 13 ) 

I F ( IS.EQ.2)  WRITE  (14)  X ,Y ,Z ,S ,P1 ,P2 ,P3 ,A1 
RETURN 


FORMAT  ( / , IX , 1 ELFAIL  l-COOR.  2-COOR.  1-STR 1 ,A3 ,3X , ' 2-STR 1 ,A 
13,'  1 2-STR 1 ,A3 ,T66 , ' PI 1 ,T77 , 'P2 1 ,T86 , ‘TAU-MAX ‘ ,T99 , 1 A1 ' ,T109 , ‘PRE 

2S' ,T120, 'E(T) ‘ ,/) 

FORMAT  ( 1111 ,5X ,20A4/5X , 1 _ ELEMENT  STRESS _ ‘,'CCW  ROTATION  IS 

1  POSITIVE _ ‘/5X,‘ _ 1-DIRECTI ON  IS  RADIAL  AND  2-DIRECTION  IS  AX 

2IAL  FOR  AN1 ,1X, ' AX  I  SYMMETRIC  PROBLEM _ '//) 

FORMAT  (/, IX, 'ELFAIL  l-COOR.  2-COOR  1-STR ', A3 ,3X , 1 2-STR ', A3 , ' 

1  3-STR ' ,A3 , 1  12-STR' ,A3,T77, 'PI ' ,T88, ' P2 ' ,T99 , ' P3 ' ,T108, ' A1 ' ,T1 

216,  'TEMP', T127 , ' E(T) ' ,/) 

FORMAT  (1X,I4,1X,A1,1P2E10.2,1P7E11.3,F7.2,1P2E11.3) 

FORMAT  (IX ,  14 ,A1 ,6E11 . 3/6X ,5E11 .3 ) 


SUBROUTINE  SETMPT  (NEXT, NOW, INCR.MTOT, TITLE) 


IMPLICIT  REAL*8( A-H ,0-Z ) 

COMMON  /SIZE/  MAX 

....SET  THE  STORAGE  POINTER  NEXT  =  NOW  +  INCR  AND  CHECK  FOR  SUFFIC 
MEMORY . 

CALLED  BY  -  ACNTRL 

MTOT=MTOT+INCR 
NEXT=NOW+INCR 
NN=NEXT-L 
N2=2*(NN/2 ) 

IF(NN .EQ.N2 )  GO  TO  1 

NEXT=NEXT+1 

MTOT=MTOT+l 

IF(MTOT.LE.MAX)  RETURN 
WRITE  (6,2)  TITLE, MTOT.MAX 
STOP 

FORMAT  (///5X,'... FATAL  ERROR. .. INSUFFICIENT  STORAGE  IN  BLANK  COMM 
ION 75X, 'REQUIRED  TO  EXECUTE  *H',A4,'*  IS '  ,18, 1  BUT  AVAILABLE  ONLY 
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SUBROUTINE  SHAPE2  (XL  , IX ,NDM,NEl ,XSJ ,M,N ,NPT) 


CALLS 


SHAP22 


SHAPE  FUNCTION  ROUTINE  FOR  2-D  ISOPARAMETRIC  ELEMENTS 
CALLED  BY  SUBROUTINE  ELEMNT 
IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON  /GQDRT/  GPT(4,4)  ,WT(4,4) ,SHP (4 ,9 ) ,LX (4 ) ,LE (4) 
DIMENSION  IX ( 1 ) ,  XL ( NDM , 1 ) ,  XS(3,2) 

X=GPT(M,NPT ) 

E=GPT(N,NPT) 

FORM  4-NODE  QUADRILATERAL  SHAPE  FUNCTION 
DO  1  1=1,4 

SHP(3,I)=(l.+LX(I)*X)*(l.+LE(I)*E)/4. 

SHP (1 , I )=LX( I )*( 1 ,+LE ( I )*E )/4. 

SHP (2 ,1 )=LE ( I )*(1 .+LX( I )*X)/4. 

IF ( NEL .GE .4 )  GO  TO  3 

FORM  TRIANGULAR  SHAPE  FUNCTIONS  BY  ADDING  3-RD  AND  4-TH 
DO  2  1=1,3 

SHP ( I ,3 )=SHP ( I ,3 )+SHP (1,4) 

ADD  QUADRATIC  TERMS  IF  NECESSARY 
I F ( MEL  - GT . 4 )  CALL  SHAP22  (X ,E ,SHP ,NEL , IX ) 

CONSTRUCT  JACOBIAN  XJ  AND  DETERMINANT  OF  XJ 
DO  4  1=1, NDM 
DO  4  J=1 ,2 
XS( I ,J )=O.ODO 
DO  4  K  =  1  ,NEL 

XS(I,J)=XS(I,J )+XL ( I ,K )*SHP (J ,K ) 

IF(NDM.NE.3)  GO  TO  5 

COMPUTE  SURFACE  AREA  AND  RETURN  AS  DETJ 

A1=XS(2,1 )*XS(3,2)-XS(3,1 )*XS(2,2) 

A2=XS(3 ,1 )*XS( 1 ,2 )-XS( 1 , 1 )*XS(3 ,2 ) 

A3=XS( 1 ,1 )*XS(2 ,2)-XS(2 ,1 )*XS(1 ,2 ) 
DETJ=DSQRT(A1*A1+A2*A2+A3*A3) 

RETURN 

DETJ=XS(1,1)*XS(2,2)-XS(1,2)*XS{2,1) 

COMPUTE  DNI/DX  AND  STORE  IN  SHP (1,1) 

COMPUTE  DNI/DY  AND  STORE  IN  SHP (2, I) 

DO  6  I =1 .NEL 

TEMP=(XS(2 ,2 )*SHP (1,1 )— XS (2,1 )*SHP(2, I ) )/DETJ 
SHP (2 , I )= ( -XS( 1 ,2 )*SHP (1 , I )+XS( 1 ,1 )*SHP (2 , I ) )/DETJ 
SHP (1,1 )=TEMP 

XSJ=DETJ*WT(M,NPT)*WT(N,NPT) 


RETURN 


o  o  o  o  o  o 
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SUBROUTINE  SHAP22  (S,T,SHP ,NEL , IX ) 


CALLED  BY  SHAPE2 

ADD  QUADRATIC  TERMS  AS  NECESSARY 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  IX ( 1 ) ,  SHP (4 ,9 ) 

C 

S2=(1.0D0-S*S)/2.0DO 
T2= (1 .ODO-T*T)/2,ODO 
DO  1  1=5, NEL 
DO  1  J  =  1 ,3 

1  SHP (J , I )=O.ODO 

C  ....  MIDSIDE  NODES  (SERENDIPITY) 

IF(IX(5). EQ . 0 )  GO  TO  2 
SHP ( 1 ,5)=-S*(1.0D0-T) 

SHP (2 ,5 )=-S2 

SHP (3 ,5 )=S2*( 1 .ODO-T ) 

2  I F ( NEL . LT . 6 )  GO  TO  8 
IF ( IX(6).EQ.O)  GO  TO  3 
SHP (1 ,6 )=T2 

SHP (2 ,6 )=-T*( 1 .ODO+S ) 

SHP (3 ,6 )=T2*( 1 .ODO+S ) 

3  IF ( NEL . LT . 7 )  GO  TO  8 
IF( IX(7 ) .EQ.O)  GO  TO  4 
SHP ( 1 ,7 )=-S*( 1 .ODO+T ) 

SHP(2,7)=S2 

SHP (3 ,7 )=S2*(1 .ODO+T ) 

4  IF(NEL .LT.8)  GO  TO  8 

I F ( I X (8 ) . EQ . 0 )  GO  TO  5 

SHP ( 1 ,8 )=-T2 

SHP (2 ,8 )=-T*( 1 .ODO-S  ) 

SHP (3 ,8)=T2*( 1 .ODO-S ) 

C  ....  INTERIOR  NODE  (LAGRANGIAN) 

5  I F ( NEL - LT . 9 )  GO  TO  8 
IF( IX (9 ) .EQ.O )  GO  TO  8 
SHP ( 1 ,9)=-S*T2 

SHP (2 ,9 )=-T*S2 
SHP (3 ,9 )=4.0D0*S2*T2 

C  ....  CORRECT  EDGE  NODES  FOR  INTERIOR  NODE  (LAGRANGIAN) 
DO  7  J=1 ,3 
DO  6  1=1,4 

6  SHP (J , I  )=SHP (J , I )-0.25D0*SHP(J ,9 ) 

DO  7  1=5,8 

7  IF(IX(I).NE.O)  SHP ( J ,1  )=SHP(J ,1 )-0.5D0*SHP(J ,9) 

C  ....  CORRECT  CORNER  NODES  FOR  PRESENCE  OF  MIDSIDE  NODES 

8  K=8 

DO  10  1=1,4 
L  =  I+4 

DO  9  J  =  1 ,3 


o  o  o 


SUBROUTINE  ZERO  (A,N) 

C  . . 

IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  A( 1 ) 

CALLED  BY  -  ACNTRL ,FORMEL 

00  1  1  =  1 ,  M 
1  A( I ) =0 . ODO 
RETURN 
END 


oooooo  oo  o  o  ooooooooo 
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SUBROUTINE  CONSOL  (X ,T, IX, IE , ID, IDL ,IXL ,KPRS,XL ,DL ,SPU ,UL ,FL ,TL , 


1W, F,R,D,AK, TIMES, NDT,N0PT,NPRINT,NEN,NEN1,  NPN , 

2NDF,NDM,MEQ, I  BAND, NUMEL.NUMNP, NSTEPS, NST.NRD.BND) 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /ANGLE/  ANG ( 4 ) ,DC(3) ,ST(2)  ,IEPS,LINE,TWR 
LOGICAL  BACK ,F0RW,BND ,SIG ,TWR 

DIMENSION  X (NDM,1 ) ,  T(l),  IX(NEN1,1),  IE (4,1 ) ,  ID(NDF.l),  IDL(l), 
1IXL(1),  KPRS(l),  XL  (1 ) ,  DL(1),  SPU(NST.l),  UL(1),  FL(1),  TL(1),  W( 
2NDF.1),  F(NOF.l),  0(1),  AK(MEQ .IBAND ) ,  R(NDF ,NUMNP ) ,  TIMES(l),  NOT 
3(1),  N0PT(1),  NPRINT(l) 

FILE  ASSIGNMENT 

15  -  ELEMENT  STIFFNESS 

16  -  GLOBAL  STIFFNESS  (MEQ, IBAND) 

17  -  GLOBAL  INITIAL  FORCE  (NDF.NUMNP) 

SIG=. FALSE. 

CORRECT  BCODE  IN  THE  8-4  ELEMENT 

IF(NPN.EQ.NEN)  GO  TO  2 

NP1=NPN+1 

DO  1  M=1 ,NUMEL 

DO  1  I=NP1 ,NEN 

J  =  I X ( I ,M) 

ID(3,J)=1 

1  CONTINUE 

2  CONTINUE 

...INPUT  TIME  FACTORS 

READ  9 ,  ( T I  ME  S ( I ) ,NDT ( I ) ,NOPT ( I ) ,NPRINT( I ) ,KPRS( I ) , 1=1 ,NSTEPS ) 
PRINT  10,  NEN ,NPN 
PRINT  11 

PRINT  12,  ( I .TIMES ( I ) ,NDT ( I ) ,NOPT ( I ) ,NPRINT( I ) ,KPRS( I ) , 1=1 .NSTEPS ) 

_ ADD  EXTERNALLY  APPLIED  CONC.  NODAL  LOADS  TO  R  AND  INITIATE  W 

DO  3  N=1 ,NUMNP 
DO  3  J=1  ,NDF 

IF(ID(J,N).NE.O)  GO  TO  3 
R(J,N)=R(J,N)+F(J,N) 

W(J  ,N)=O.ODO 

3  CONTINUE 

....STORE  THE  GLOBAL  STIFFNESS  MATRIX  TO  BE  MULTIPED  BY 
ALPHA*DT  IN  SUBROUTINE  TMSTP  AMD  THE  INITIAL  LOAD  VECTOR 
WRITE  (16)  AK,R 


....CONSOLIDATION  COMPUTATIONS 
DO  8  LT=1, NSTEPS 
BACK=. FALSE. 


r. 


LT 1 =LT — 1 
ST E  P  =  T I  ME  S  ( L  T ) 

IF(LT.GT.l)  STEP =T I MES ( LT )-TIMES(LTl ) 

DT=STEP/NDT (LT ) 

I F ( NOPT ( LT ) .GT.4)  ALPHA=1.0D0+(1.0D0/DT)-(1.0D0/(DL0G(1.0D0+DT) ) ) 

IF(NOPT(LT).EQ.l)  ALPHA=0.5D0 

IF ( NOPT ( LT  )  .EQ.2 )  ALPHA=0. 666600 

IF(NOPKLT) .EQ.3)  ALPHA=0 .87100 

I F ( MOPT ( LT ) .EQ.4)  ALPHA=01.D0 

IF(LT.EQ.l)  GO  TO  4 

CONTINUE 

MULTIPLY  K2  IN  AK  BY  ALPHA  *  DT ,  APPLY  GEOMETRIC  BOUNDARY 
CONDITIONS,  AND  SAVE  R  (LOAD  VECTOR) 

CALL  TMSTP  ( ALPHA ,DT , R ,W ,F ,AK , IBAND ,MEQ .NUMNP ,NUMEL , IX , ID ,NEN1 ,NDF 
1  ,SPU ,NST ,NPN  ) 

....STORE  THE  GEOMETRICALLY  MODEFIED  GLOBAL  LOAD 
VECTOR  TO  BE  UGDATED  IN  SUBROUTINE  TMINC 
REWIND  17 
WRITE  (17)  R 

TRAINGULIZE  THE  STIFFNESS  MATRIX 

CALL  BAMSOL  ( AK ,R,R ,MEQ , IBAND , BACK ,FORW ) 

UPDATE  THE  LOAD  VECTOR  R,  FOR  EVERY  TIME  INCREMENT 
NDTS=NDT ( LT ) 

SET  AN  OUTPUT  COUNTER 
NO=MPRINT(LT ) 

DO  7  LDT=1,NDTS 

CALL  TMINC  (R.W.SPU , ID , IX ,NDF,NEN ,NEN1 , NUMNP ,NUMEL ,NST ,NPN, ALPHA, D 
IT) 

SOLVE  FOR  DISPL.  AND  PORE  PRESS.  CORRESPONDING  TO  R 
BACK=.TRUE. 

CALL  BANSOL  ( AK , R ,R ,MEQ , IBAND .BACK ,FORW ) 

DO  5  1  =  1,  NUMNP 
DO  5  J=1,NDF 
W(J , I  )  =  R(J  ,  I  ) 

TIME=TIMES(LT )+LDT*DT-NDTS*DT 
IF (LDT.EO.NO. OR. LDT.EQ.NDTS )  GO  TO  6 
GO  TO  7 

CALL  PRTDIS  ( ID ,X ,R,F ,NDM,NDF , .FALSE . ,TWR ..TRUE . ,LT ,LDT .TIME ,BND ) 
NO=NPRINT(LT )+N0 
I F ( KPRS ( LT ) .EQ.O )  GO  TO  7 
LINE  =0 

CALL  FORMEL  (X ,XL ,D ,DL  ,SPU  ,UL  ,FL  ,TL  ,R ,F ,T , IX , IE , ID, IDL , IXL , IXL ,NDF 
1  ,NDM,NEN1 ,NST,5 ,AK ,AK ,R,SIG ,SIG ,SIG ,SIG ,SIG ,NRD,NOW, IBAND ,BND ) 
CONTINUE 
CONTINUE 
RETURN 


FORMAT  (1F10. 0,415) 

FORMAT  ( 1H1/15X , 'CONSOLIDATION  ANALYSIS' //10X, 'NUMBER  OF  ELEMENT  N 
10DES  =' ,I5/10X, 'NUMBER  OF  PRESSURE  NODES  =  * ,15//) 

FORMAT  (6X, 'STEP  1 ,7X ,' ELLAPSED ' ,2X , 'NUMBER  OF ' ,6X, 'TIME ' ,4X, 'OUTPU 
IT ',3X,' STRESS '/4X, 'NUMBER ',11X, 'TIME' ,4X, 'TIM  INC ' ,2X ,' INTERPOL ' ,3 
2X,' COUNTER* ,3X, 'OUTPUT' ) 

FORMAT  (110,1015.5,4110) 


SUBROUTINE  TMSTP  (ALPHA, DT,R,W,F,AK, IBAND, 

1MEQ,NUMNP ,NUMEL , IX , ID ,NEN1 ,NDF,SPU ,NST ,NPN ) 

IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  R(NDF.NUMNP),  W(NDF.l),  F(NDF,1 ) ,  AK (MEQ , IBAND ) ,  ID (NDF , 
11),  IX (NEN1 ,1 ) ,  CPP (8,8),  SPU(NST.NST) 

REWIND  16 
READ  (16)  AK,R 

MULTIPLY  K2  IN  THE  GLOBAL  STIFFNESS  BY  ALPHA*DT 

DO  2  1=1 ,NUMNP 

I I=NDF*I 

DO  1  J=1 ,NUMNP 

JJ=NDF*J 

KK=JJ-I 1+1 

IF(KK.LE.O. OR. KK.GT. IBAND)  GO  TO  1 
AK ( 1 1 ,KK )=AK (II ,KK )*ALPHA*DT 
CONTINUE 
CONTINUE 

ADD  CPP  TO  AK  GLOBALLY  AND  TO  SPU  LOCALLY 
REWIND  15 
DO  4  M=1 ,NUMEL 

READ  ( 15 )  NSL , ( (SPU ( I ,J ) , J  =  1 ,NSL ) , I  =  1 ,NPN ) , ( (CPP ( I ,J ) ,J=1 ,NPN ) , 1=1 
1  ,NPN ) 

DO  3  1  =  1, NPN 
I  I=NDF*IX ( I ,M ) 

DO  3  J=1  ,NPN 
JJ=NDF*IX(J,M) 

KK=JJ-II+1 

IF(KK.LE.O)  GO  TO  3 

AK  ( 1 1  ,KK  ) = AK  (II  ,KK  )-CPP  ( I , J  ) 

CONTINUE 

CONTINUE 

INTRODUCE  KINEMATIC  CONSTRAINTS 
DO  6  M=1 ,NUMNP 
DO  5  J  =  1  ,NDF 

IF( ID (J ,M) .EQ.O)  GO  TO  5 
IDF=(M-1 )*NDF+J 

CALL  GEOMBC  (F(J ,M) , IDF, IBAND, MEQ, R,AK ) 

CONTINUE 

CONTINUE 


SUBROUTINE  TMINC  (R,W,SPU , ID , IX ,NDF ,NEN ,NEN1 .NUMNP ,NUMEL ,NST ,NPN 


1ALPHA.DT) 

IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  R(NDF, NUMNP),  W(l),  ID(NDF.l),  IX(NENl.l),  SPU (NST.NST ) , 
1  CPP (8 ,8 ) 

REWIND  15 
REWIND  17 
READ  (17)  R 
DO  4  M=1 ,NUMEL 

READ  (15)  NSL , ( ( SPU ( I ,J ) ,J  =  1 ,NSL ) , 1=1 ,NPN ) ,  ( (CPP ( I  ,J )  ,J  =  1 ,NPN ) ,  1  =  1 
1 ,NPN ) 

NEL=NSL/NDF 
DO  3  1=1, MPN 
J  =  I X ( I ,M) 

IF(ID(3,J).NE.O)  GO  TO  3 
DO  2  K=1 ,NEL 
IF(K.GT.NPN)  GO  TO  1 
KK=K*NDF 

SPU ( I ,KK )=SPU ( I ,KK )*( (ALPHA-1 .DQ)*DT )-CPP ( I ,K ) 

CONTINUE 

N=IX(K,M) 

DO  2  L=1 ,NDF 
KLOC=(K-l )*NDF+L 
KGLB= (M-l )*NDF+L 

R(3 ,J )=R(3,J )+SPU( I ,KLOC )*W(KGL8) 

CONTINUE 

CONTINUE 

RETURN 


SUBROUTINE  ELMTOl  ( IX ,XL  ,D  ,FL  ,S,UL ,TL .NDM.NDF ,NST ,NEL , ISW 


1IL,LN,NG) 

CALLED  BY  ELMLIB 

CALLS  SHAPE2 .PVALUE ,GAB05I 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /CDATA/  0 , HEAD (20 ) ,NUMNP ,NUMEL,NUMMAT ,NEN ,NEO ,MEQ .NTAPE ,NPN 
COMMON  /ANGLE/  ANG (4 )  ,DC{3 )  ,ST(2 ) , IEPS .LINE ,EWR 
COMMON  /GQDRT /  GPT{4,4) ,WT(4,4) ,SHP(4,9) ,LX(4) ,LE (4) 

COMMON  /STRES/  EPS(6  )  ,SIG(6 ) 

DIMENSION  D(l),  XL(NDM.l),  FL(NDF,1),  S(NST.NST),  UL(NDF,1),  TL ( 1 ) 
1,  IX ( 1 ) ,  SS(3,8),  CPP (8,8) 

NGP=NG 

NDU=NDF-1 

GO  TO  (1,1,15,16,21),  ISW 
CONTINUE 


....COMPUTE  ELEMENT  STIFFNESS  MATRIX  BY (NGP*NGP ) 

NSL=NEL*NDF 

DO  2  1=1, NPN 

DO  2  J  =  1  ,NPN 

CPP ( I  ,J )=0.0D0 

IF(MEL.EQ.8)  GO  TO  3 

CALL  GABOSI  ( IX ,XL ,D ,S ,CPP , NDM.NDF ,NST ,6 ,NPN ,NGP ) 

GO  TO  14 

AK1=-D(24) 

AK2=-D (25 ) 

C0MP=D(22) 

DO  10  M=1 ,NGP 
DO  10  N  =  1  ,NGP 

CALL  SHAPE2  (XL , IX ,NDM,NEL ,XDJ ,M,N ,NGP ) 

DO  4  1=1, NDU 
DO  4  J  =  1  ,NEL 
SS ( I , J )=SHP ( I ,J ) 

CALL  SHAPE2  (XL,IX,NDM,NPN,XPJ ,M,N,NGP ) 

....COMPUTE  STIFFNESS  MATRIX  OF  THE  SOLID  PHASE 

COLUMN 

Jl  =  l 

DO  6  J=1,NEL 
W1J=SS{ 1 ,J )*XDJ 
W2J=SS( 2 , J )*XDJ 
K 1  =J  1 
ROW 

DO  5  K=J  ,NEL 

S(J1,K1)=S(J1,K1 )+WlJ*SS( 1 ,K ) 
S(J1,K1+1)=S(J1,K1+1)+W1J*SS(2,K) 
S(J1+1,K1)=S(J1+1,K1 )+W2J*SS( 1 ,K ) 
S(J1+1,K1+1)=S(J1+1,K1+1)+W2J*SS(2,K) 


o  o 
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5  K1=K1+NDF 

6  J1=J1+NDF 

COMPUTE  STIFFNESS  MATRIX  OF  THE  FLUID  PHASE 
DO  7  1=1  ,NPN 
1 1  =  1 *NDF 

WI I =SHP (3,1) *XPJ  *COMP 
W1I=SHP (1,1 )*XPJ*AK1 
W2I=SHP (2,1 )*XPJ*AK2 
DO  7  J  =  I  ,NPN 
JJ=J*NDF 

CPP ( I ,J )=CPP ( I ,J )+WI I*SHP (3  ,J ) 

7  S ( 1 1 ,JJ )  =  S ( 1 1 , J J )+Ml I *SHP ( 1 , J )+W2I *SHP (2 ,J ) 

C 

C  COMPUTE  COUPLING  STIFFNESS  MATRIX 

11=1 

DO  9  1=1, NEL 
Wl I=SS( 1 , I )*XDJ 
W2I=SS(2,I )*XDJ 
00  8  J  =  1  ,NPN 
JJ=J *NDF 

S( 1 1 ,JJ )  =  S( 1 1 ,JJ )+Wl I*SHP (3,J ) 

8  S( 1 1+1 ,JJ )  =  S( 1 1+1 ,JJ )+U2I*SHP {3 ,J ) 

9  1 1  =  I I+NDF 

10  CONTINUE 

DO  11  J=1 ,NSL ,NDF 
DO  11  K=J ,NSL ,NDF 
W11=S( J ,K ) 

W12=S( J ,K+1 ) 

U21=S(J+1,K) 

W22=S( J+l ,K+1 ) 

S(J ,K )=W11*D(1 )+W22*D(4)+W12*D(5 )+W21*D(5 ) 

S( J ,K+1 )=W12*D(2 )+W21*D(4 )+Wll*D(5 )+W22*D(6 ) 

S(J+1 ,K )=W21*D(2)+W12*D(4)+W22*D(6)+W11*D(5) 
S(J+1.K+1)=W22*D(3)+W11*D(4)+W21*D(6)+W12*D(6) 

S(K,J)  =  S(J,K) 

S(K,J+1)=S(J+1,K) 

S(K+1 ,J )=S(J ,K+1 ) 

11  S(K+1,J+1)=S(J+1,K+1) 

11  =  1 

DO  13  1=1, NEL 
00  12  J  =  1  ,NPN 
JJ=J *NDF 

S(JJ , 1 1 )=S( 1 1  ,JJ  ) 

S{JJ ,11+1 )=S{ 1 1  +  1  ,JJ  ) 

DO  12  K  =  1  ,NPN 
KK=K*NDF 

CPP  (K  ,J  )=CPP  (J  ,K  ) 

12  S(KK,JJ)=S(JJ,KK) 

13  1 1  =  1 I+NDF 

14  NPL=NDF*NPN 

WRITE  (15)  NSL , ( (S( I  ,J  )  ,J  =  1  NSL ) , I=NDF ,NPL ,NDF ) ,  ( (CPP ( I  ,J  )  ,J=1  ,NPN 


1 ) , I =1 ,NPM  ) 

RETURN 

CONTINUE 

....GET  FLUID  AND  SOLID  80DY  FORCES 
AK  2 = D  ( 2  5 ) 

R0WF=-D(26 ) 

DO  20  M=1 ,NGP 
DO  20  N=1 ,NGP 

CALL  SHAPE2  (XL , IX .NDM.NEL ,XDJ ,M,N ,NGP  ) 
DO  17  1=1, NEL 
SS(3,I )=SHP(3,I ) 

CALL  SHAPE2  (XL , IX .NDM.NPN ,XPJ ,M,N ,NGP ) 
D2=XDJ*D(23) 

DO  19  1  =  1, NEL 
D1=D2*SS(3,I ) 

DO  18  J  =  1  ,NDU 

FL(J,I)=FL(J,I)+D1*DC(J) 

IF(I.GT.NPN)  GO  TO  19 

FL(3,I)=FL(3,I )+AK2*R0WF*SHP (2 , I )*XPJ 

CONTINUE 

CONTINUE 

RETURN 

CONTINUE 


- COMPUTE  STRESSES  AT  NGP  GAUSS  PTS. 

Z=0 .000 
DO  24  M=1 ,NGP 
DO  24  M=1,NGP 

CALL  SHAPE2  (XL , IX , NDM.NEL ,XSJ ,M,N , NGP ) 

DO  22  1=1,3 

EPS ( I  )=0 .ODO 

T=O.ODO 

X=O.ODO 

Y=O.ODO 

DO  23  J=1,NEL 

T=T+SHP(3,J)*TL(J) 

X=X+SHP (3 ,J )*XL (1 ,J ) 

Y=Y+SHP ( 3 , J  )*XL ( 2  ,J  ) 

EPS(l)  =  STRAIN  IN  X-DIRECTION  ,  EPSILON-XX 

EPS(2 )  =  STRAIN  IN  XY-DIRECTION  ,  GAMMA-XY 

EPS ( 3  )  =  STRAIN  IN  Y -DIRECTION  ,  EPSILON-YY 

EPS ( 1  )=EPS(1 )+SHP (1 ,J )*UL( 1 ,J ) 
EPS(3)=EPS(3)+SHP(2,J)*UL(2,J) 

EPS(2 )=EPS(2 )+SHP (2 ,J )*UL ( 1 ,J }+SHP (1 , J )*UL (2 ,J ) 
EPS { 1 ) =EPS ( 1 )-D(8)*T 
EPS(3 )=EPS(3 )-D (9 )*T 


SIG(l)  =  STRESS  IN  X-DIRECTION 
S I G ( 2 )  =  STRESS  IN  XY-DIRECTION 
S I G ( 3  )  =  STRESS  IN  Y-DIRECTION 


SIGMA-XX 

TAU-XY 

SIGMA-YY 


CALCULATE  EFFECTIVE  STRESSES 


SIG(1)=D(1)*EPS(1)+D(2)*EPS(3) 

S I G ( 2 )=D(4)*EPS(2 ) 

SIG(3 )=Q(2 )*£PS(1 )+D(3 ) *EPS(3 ) 

CALL  PVALUE  (SIG,1,2,X,Y,Z,LN,IL,TE,D) 

....CALCULATE  TOTAL  STRESSES 

CALL  SHAPE2 (XL , IX ,NDM ,NPN ,XSJ ,M,N ,NGP ) 

DO  49  J=1 ,NPN 

SIG(1)=SIG(1 )+SHP (3,J )*UL (3, J ) 

49  SIG(3)=SIG(3)+SHP (3,J)*UL(3,J) 

CALL  PVALUE (SIG, 1, 2, X,Y,Z,LN,IL,TE,D) 

CONTINUE 

RETURN 

END 
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SUBROUTINE  GEOMBC  (U  ,N  ,  I8AND ,MEQ,R ,AK ) 

....THIS  SUBROUTINE  APPLIES  THE  KINEMATIC  CONSTRAINTS  FOR  EVERY  D 
....CALLED  BY  TMSTP 
IMPLICIT  REAL*8(A-H ,0-Z ) 

DIMENSION  AK ( ME 0,1) ,  R(l) 

DO  2  M=2 , IBAND 
K=N-M+1 

IF(K.LE.O)  GO  TO  1 
R  ( K  )  =  R  ( K  )  -  AK  ( K  ,  M )  *U 
AK (K ,M)=O.ODO 
K=N+M-1 

IF(K.GT.MEQ)  GO  TO  2 

R(K  )  =  R  (K  )-AK  (N  ,M )  *U 

AK(N,M)=O.ODO 

CONTINUE 

AK  (N  ,1 )  =  1 .ODO 

R(M  )=U 

RETURN 
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SUBROUTINE  GABOSI  ( I X , XL ,D,S,CPP,NDM,NDF,NST,NEL,NPN,NGP ) 


CALLS  SHAPEX 

IMPLICIT  REAL*8(A-H ,0-Z ) 

COMMON  /GQDRT/  GPT(4,4 ) ,WT(4 ,4) ,SHP (4 ,9) ,LX(4) ,LE (4) 
DIMENSION  D(l),  XL(NDM,1),  S(NST.NST),  IX(1),  CPP(8,8) 
1SUU (12,12),  SPU(12,4) ,  LP (4) ,  L0(8) 

_ _ INITIALIZE  STIFFNESS  MATRICES 

NDU=NDF-1 
NPL=NPN*NDU 
NSL»NEL*NDU 
DO  3  1=1, NSL 
DO  1  J=1 ,NSL 
SUU( I ,J ) =0 . ODO 
DO  2  K=1  ,NPN 
SPU ( I ,K )=0.0D0 
CONTINUE 
DO  4  1  =  1  ,NPN 
DO  4  J=1  ,NPN 
SPP ( I  ,J  )=0.0D0 

....COMPUTE  ELEMENT  STIFFNESS  MATRIX  BY (NGP*NGP ) 
AK1=-D{24) 

AK2=-D(25) 

C0MP=D(22) 

DO  10  M=1 ,NGP 
DO  10  N=1 ,NGP 

CALL  SHAPEX  (XL , IX .NDM.NEL ,XDJ ,M,N ,NGP ) 

....COMPUTE  STIFFNESS  MATRIX  OF  THE  SOLID  PHASE 

COLUMN 

Jl  =  l 

DO  6  J=1  ,NEL 
W1J=SHP ( 1 ,J )*XDJ 
W2J=SHP (2  ,J )*XDJ 
K 1  =J  1 
ROW 

DO  5  K=J  ,NEL 

SUU ( J 1 ,K1 )=SUU(J1 ,K1 )+WlJ*SHP(l ,K ) 

SUU (J1 ,K1+1 )=SUU (J1 ,K1+1 )+WlJ*SHP (2  ,K ) 

SUtKJl+l ,K1 )=SUU (J 1+1 ,K1 )+W2J*SHP (1 ,K ) 

SUU (J 1+1 ,K1+1 )=SUU ( J 1+1 ,K1+1 )+W2J*SHP(2,K ) 

K 1 =K 1 +NDU 

J1=J1+NDU 

COMPUTE  STIFFNESS  MATRIX  OF  THE  FLUID  PHASE 

DO  7  1  =  1, NPN 

WII=C0MP*SHP(3,I)*XDJ 

W1I=SHP(1 ,1 )*XDJ*AK1 

W2I=SHP(2,I )*XDJ*AK2 


,  SPP (4 ,4) , 


00  7  J=I,NPN 

CPP ( I ,J )=CPP ( I ,J )+WI I *SHP (3 ,J ) 

7  SPP  ( I , J  )  =  SPP  ( I , J  )+Wl  I  *SHP  ( 1 ,  J  )+W2I *SHP (2  ,J  ) 

C 

C  COMPUTE  COUPLING  STIFFNESS  MATRIX 

11  =  1 

DO  9  1=1, MEL 
W1 I =SHP (1,1) *XDJ 
W2I=SHP(2,I )*XDJ 
00  8  J  =  1  ,NPN 

SPU  ( 1 1 , J )  =  SPU { 1 1 , J )+Wl I *SHP ( 3  ,J  ) 

8  SPU (II+1,J)  =  SPU(II  +  1,J )+W2I*SHP(3,J  ) 

9  1 1  =  1 1+NDU 

10  CONTINUE 

DO  11  J  =  1  ,NSL ,NDU 
DO  11  K=J ,NSL  ,NDU 
W11=SIJU  ( J  ,K  ) 

W12=SUU ( J  ,K+1 ) 

W21=SUU(J+1  ,K ) 

W22=SUU (J+l  ) 

SUU ( J ,K)=Wl 1*0(1 )+W22*D(4)+W12*D(5)+W21*D(5) 
SUU(J,K+l)=W12*D(2)+W21*D(4)+Wll*D(5)+W22*D(6) 
Sill) (J+l  ,K)=W21*D(2 )+W12*D(4)+W22*D(6)+Wll*D(5) 
SUU (J+l ,K+1 )=W22*D(3)+W11*D(4)+W21*D(6)+W12*D(6) 
SUU(K,J  )=SUU(J,K) 

SUU (K , J+l )=SUU (J+l ,K ) 

SUU (K+l , J )=SUU ( J ,K+1 ) 

11  SUU(K+1  ,J+1  )=SUIJ (J  +  l  ,K+1 ) 

DO  12  1=1, NPN 

DO  12  J  =  1  ,NPN 
SPP (J  , I  )=SPP(I  ,J) 

12  CPP(J,I)=CPP(I,J) 

C  ....STATIC  CONDENSATION 
C 

IF(MPN.EQ.NEL)  GO  TO  18 
INTER=NDU*(NEL-NPN ) 

FAC=1 .DO 

IF (COMP .NE .0.000)  FAC=2.D0 
DO  17  N=l, INTER 
L=MSL-N 
M=L  +  1 

PI VOT=SUU (M ,M) 

DO  13  1=1, NPN 
FP=FAC*SPU(M,I ) /PIVOT 
DO  13  J  =  1  ,NPN 

13  CPP ( I , J )=CPP ( I ,J )+FP*SPU (M ,J ) 

DO  16  J=1,L 

FU=SUl)  (M,J ) /PIVOT 
DO  14  K  =  1  ,NPN 

14  SPU ( J  ,K  )=SPU(  J  ,K  )-Fll*SPll (M,K ) 

DO  15  I =J  ,L 

SIJU  ( I ,  J  )  =SUU  ( I ,  J  )-FU*SUU  ( I  ,M ) 


15  SUU ( J , I )  =  SUU ( I  »J  ) 

16  CONTINUE 

17  CONTINUE 

18  CONTINUE 

....RELOCATE  ELEMENT  MATRICES 
MM=0 

DO  19  I =NDU  ,NPL  ,NDU 
LD ( I  — 1 )=MM+1 
LD ( I )=MM+2 

19  MM=MM+NDF 

DO  20  1=1  ,NPN 

20  LP ( I )=NDF*I 
00  21  1=1, NPL 
LI=LD( I ) 

DO  21  J  =  1  ,NPL 
LJ=LD( J ) 

21  S ( L I ,LJ )=SUU ( I , J  ) 

DO  22  1=1, NPN 

L I  — LP ( I ) 

DO  22  J=1,NPM 
LJ=LP (J ) 

22  S ( L I ,LJ )=SPP ( I , J  ) 

DO  23  1=1, NPL 
LI=LD( I ) 

DO  23  J=1,NPN 
LJ=LP  (J ) 

S(L I ,LJ )=SPU ( I  ,J  ) 

23  S(LJ,LI)=SPU(I,J) 

RETURN 
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SUBROUTINE  SHAPEX  (XL ,IX ,NDM,NEL ,XSJ ,M,M ,NPT) 


CALLED  BY  GBABOSI 

SHAPE  FUNCTION  ROUTINE  FOR  2-D  ISOPARAMETRIC  ELEMENTS 
CALLED  BY  SUBROUTINE  ELEMNT 
IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GQDRT/  GPT(4,4 ) ,WT(4,4 )  ,SHP (4,9 ) ,LX (4 ) ,LE (4 ) 
DIMENSION  IX ( 1 ) ,  XL(NDM,1 ) ,  XS(3,2) 

X=GPT(M,NPT ) 

E=GPT(N ,NPT ) 

FORM  4-NODE  QUADRILATERAL  SHAPE  FUNCTION 
DO  1  1=1,4 

SHP(3,I)=(1.+LX(I )*X)*(l.+LE(I)*E)/4. 

SHP (1,1 )=LX( I )*(1 .+LE( I )*E )/4. 

SHP ( 2 ,1 )=LE( I )*(1 .+LX( I )*X)/4. 

ADD  GHABOUSSI  INCOMPATIBLE  SHAPE  FUNCTIONS 

SHP(3,5)=i.DO-X*X 

SHP ( 1 ,5 )=-2.DO*X 

SHP(2,5)=O.ODO 

SHP (3 ,6 )=1 ,DO-E*E 

SHP ( 1 ,6 )=O.ODO 

SHP(2,6)=-2.DO*E 

CONSTRUCT  JACOBIAN  XJ  AND  DETERMINANT  OF  XJ 

DO  2  1=1, NDM 

DO  2  J=1 ,2 

XS(I,J)=O.ODO 

DO  2  K  =  1 ,4 

XS(I,J)=XS(I,J)+XL(I,K)*SHP(J,K) 


DETJ=XS (1,1) *XS (2,2)-XS(I,2) *XS (2.1) 

COMPUTE  DNI/DX  AND  STORE  IN  SHP (1,1) 

COMPUTE  DNI/DY  AND  STORE  IN  SHP (2, 1) 

DO  3  1  =  1, MEL 

TEMP=(XS(2 ,2 )*SHP (1,I)-XS(2,1 )*SHP (2 , I ) )/DETJ 
SHP (2 , I )= (-XS(1 ,2 )*SHP (1,1 )+XS(l ,1 )*SHP (2, I ) )/DETJ 
SHP (1,1 )  =  TEMP 

XSJ=DETJ*WT(M,NPT )*WT(N ,NPT ) 

RETURN 

END 


ooo  oooooo 


98 


SUBROUTINE  DYNAMC  (X ,T , IX , IE , ID, IDL , IXL ,KPRS ,XL ,DL ,SK ,UL ,FL ,TL , 

C  . . . . . . . . 

1U0 , F, R,D,AK, TIMES, NOT .NPRINT, THETA, GAMA, BETA, VO, AO, Ul, VI, A1 , 

2NEN ,NEN1 ,NDF ,N0M ,MEQ , IBAND .NUMEL ,NUMNP .NSTEPS ,NST ,NRD,BND ) 

C 

IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON  /ANGLE/  ANG ( 4 ) ,DC(3) ,ST{2) , IEPS, LINE ,TWR 
LOGICAL  BACK ,FORW,BND,SIG,TWR 

DIMENSION  X (NDM, I ) ,  T(l),  IX(NEN1,1),  IE (4,1 ) ,  ID(NDF,NUMNP ) ,  IDL ( 
11),  IXL ( 1 ) ,  KPRS(l),  XL ( 1 ) ,  DL  ( 1 ) ,  SK(NST.l),  UL(1),  FL ( 1 ) ,  TL ( 1 ) , 

2  UO(NDF.l),  F(NDF.l),  D(l),  AK (MEQ, IBAND ) ,  R(NDF.NUMNP) ,  TIMES! 1 ) , 

3  NDT(l),  NPRINT ( 1 ) ,  THETA ( 1 ) ,  GAMA ( 1 ) ,  BETA ( 1 ) ,  VO(NDF.l),  AO ( NDF , 
41),  Ul (NDF,1 ) ,  VI (NDF, 1 ) ,  A1(NDF,1),  COEF(IO),  TRECD(20),  ACEL (20 ) 
5,  I ACEL (4 ) 

FILE  ASSIGNMENT 

15  -  ELEMENT  STIFFNESS 

16  -  GLOBAL  STIFFNESS  (MEQ, IBAND) 

17  -  GLOBAL  INITIAL  FORCE  (NDF.NUMNP) 


SIG=. FALSE. 


INPUT  INITIAL  CONDITIONS 

READ  9,  USO ,VSO ,UFO ,VFO 

PRINT  10,  USO ,VSO ,UFO ,VFO 

NDF2=NDF/2 

DO  2  1=1 ,NUMNP 

DO  1  K=1 ,NDF2 

L=NDF-K+1 

UO (K , I )=USO 

UO ( L , I )=UFO 

VO(K ,1 )=VSO 

1  VO ( L , I )=VFO 
DO  2  J  =  1  ,NDF 

IF( ID(J,I ).EQ.O)  GO  TO  2 
U0( J ,1  )=F( J  ,1 ) 

VO(J ,1 )=0.0D0 

2  CONTINUE 

READ  11,  LOAD, OMEGA, AZERO ,BZERO ,CZERO 
PRINT  12,  LOAD, OMEGA, AZERO, BZERO.CZERO 
IF(LOAD.EQ.O)  GO  TO  3 
READ  13,  LACEL , ( I ACEL ( I ) , 1  =  1 ,NOF ) 

PRINT  14,  LACEL 
PRINT  15 

PRINT  16,  ( I ACEL ( I ) , 1  =  1 ,NDF  ) 

IF(LACEL.EQ.O)  GO  TO  3 

READ  17,  NTAPE ,HDTACL 

PRINT  18,  NTAPE, NDTACL 

READ  (NTAPE)  (TRECD( I ) ,ACEL ( I ) , 1=1 .NDTACL ) 

PRINT  19,  (ACEL ( I ) ,1=1 .NDTACL ) 

3  CONTINUE 


ooo  o  o  c~«  onoo  ooo  ooo 
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C 

C  ...INPUT  TIME  FACTORS 

READ  20 .  (TIMES! I ) , THETA ( I ) ,GAMA( I ) ,BETA( I ) ,NDT ( I ) ,NPRINT( I ) ,KPRS( 
1 1 )  ,1  =  1  ,NSTEPS) 

PRINT  21 

PRINT  22,  ( I .TIMES ( I ) , THETA( I ) ,GAMA( I ) ,BETA( I ) ,NDT (I),NPRINT(I),KP 
1RS( I )  ,1  =  1 .NSTEPS) 

....STORE  THE  GLOBAL  STIFFNESS  MATRIX  TO  BE  UPDATED 
EVERY  STEP  AND  THE  INITIAL  LOAD  VECTOR 
WRITE  (16)  AK,R 

DETERMINE  THE  INITIAL  ACCELERATION  FROM  THE  EQUATION  OF  MOTION 

CALL  INACL  ( AK ,U0 ,V0 ,A0 ,R ,F , SK , IX , I D ,NEN ,NEN1 ,NST ,MEQ , I  BAND ,NUMNP , 
1NUMEL ,NDF ,BZER0 ) 

CALL  ZERO ( AO ,NDF*NUMNP ) 


....DYNAMIC  COMPUTATIONS 
DO  8  LT=1, NSTEPS 
NO=NPRINT(LT ) 

NDTS=NDT ( LT ) 

LT1=LT-1 
STEP-TIMES ( LT ) 

IF(LT.GT.l)  STEP=TIMES(LT ) -T I  ME  S ( LT 1 ) 

DT-STEP/NDT ( LT ) 

CALCULATE  THE  DYNAMIC  STIFFNESS  MATRIX,  APPLY  GEOMETRIC  BOUNDARY 
CONDITIONS,  AND  SAVE  R  (LOAD  VECTOR) 

CALL  STEPS  ( THETA (LT ) ,GAMA(LT ) ,BETA(LT ) ,DT ,R,F,AK ,SK ,ID ,COEF ,MEQ , I 
1BAND ,NUMEL ,MUMNP ,NDF ,NST ) 

....STORE  THE  GEOMETRICALLY  MODEFIED  GLOBAL  LOAD 
VECTOR  TO  BE  UBDATED  IN  SUBROUTINE  INCRS 
REWIND  17 
WRITE  (17)  R 

C 

C  TRAINGULIZE  THE  STIFFNESS  MATRIX 
BACK=. FALSE. 

FORW= .TRUE . 

CALL  BANSOL  ( AK , R , R ,MEQ , IBAND .BACK ,FORW ) 

C 

C  UPDATE  THE  LOAD  VECTOR  R,  FOR  EVERY  TIME  INCREMENT 
C 

DO  7  L DT  =  1 ,NDTS 
TIME=T1MES(LT )+LDT*DT-NDTS*DT 

CALL  INCRS  ( THETA( LT ), GAMA ( LT), BETA (LT),DT,R,UO, VO, AO, Ul.Vl.SK, ID, 
1 1 X, NDF,NUMNP,NUMEL,NST,NEN,NEN1,F,A1, OMEGA, AZERO ,BZERO ,CZERO , TIME, 
2L0AD .LACEL , IACEL ,TRECD,ACEL .NDTACL  ) 

C 

C  SOLVE  FOR  DISPLACEMENTS(R)  AT  TO+THETA*DT 


ooooo  o o  ooo 
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BACK = .TRUE . 

CALL  RANSOL  ( AK ,R , R.MEQ , IBAND .BACK ,F0RW ) 

C 

C  OBTAIN  DISPLACEMENT,  VEL .  AND  ACC.  AT  THE  END  OF  AN  INCREMENT 
C  (IE.  AT  TO+DT) 

C 

DO  4  I =1 ,NUMMP 
DO  4  J=1 ,NDF 

Ul(J,I)=R(J,I )*C0EF ( 1 )+U0 ( J , I )*C0EF(2  )+V0(J , I )*C0EF (3 )+A0 (J , I )*C0E 
IF  ( 4 ) 

VI (J, I )=(R(J,I )-U0(J,I))*C0EF{5)+V0(J,I)*C0EF(6)+A0(J,I)*C0EF(7) 

4  Al(J,I)=(R(J,I ) — UO ( J , I ))*C0EF(8)-V0(J,I )*C0EF(9)+A0(J,I )*C0EF(10) 

OBTAIN  INITIAL  CONDITION  FOR  NEXT  INCREMENT 

DO  5  1=1 ,NUMNP 
00  5  J=1 ,NDF 
U0(J,I)=U1(J,I) 

V0(J,I)=V1(J,I) 

AO ( J , I )=A1 (J , I ) 

5  CONTINUE 

CALL  ZERO(AO,NDF*NUMNP) 

CALL  ZERO(VO  ,NDF*NUMNP ) 

IF(LDT.EQ.NO.OR.LOT.EQ.NDTS)  GO  TO  6 
GO  TO  7 

6  CALL  PRTDIS  ( ID ,X ,U1 ,F ,MDM,NDF , .FALSE . ,TWR , .TRUE . ,LT ,LDT .TIME ,BND ) 

*IN  CASE  THE  VELOCITY  (VI)  IS  REQUIRED  TO  BE  PRINTED, 

^REMOVE  THE  COMMECENT  (C)  ON  THE  FOLLOWING  TWO  CARDS 

CALL  PRTDIS ( ID, X ,V1 ,F,NUM,NDF , . FALSE., TWR, .TRUE . ,LT , LOT, TIME, 

1  BND ) 

NO=NPRINT(LT )+N0 

I F ( KPRS ( LT ) . EQ . 0 )  GO  TO  7 

LINE=0 

CALL  FORMEL  (X ,XL ,D ,DL ,SPU ,UL ,FL ,TL ,U1 ,F ,T , IX , IE , ID , IDL ,  IXL  ,IXL  ,ND 
IF ,NDM,NEN1 ,NST ,5 ,AK ,AK ,R,SIG ,SIG ,SIG ,SIG ,SIG ,NRD , NOW, IBAND , BND ) 

7  CONTINUE 

8  CONTINUE 
RETURN 

C 

9  FORMAT  (4F10.0) 

10  FORMAT  (1H1/10X, ' INITIAL  CONDITIONS..  7/10X, 'SOLID  INITIAL  DISPLA 

1CEMENTS  =  '.E10.3/10X, ‘SOLID  INITIAL  VELOCITY  =  ’  .E10.3/10X , ' 
2FLUID  INITIAL  DISPLACEMENT  =  ‘ ,E10 .3/10X , ' FLUID  INITIAL  VELOCITY 
3  =  1 ,E10.3///) 

11  FORMAT  (I5.4F10.0) 

12  FORMAT  (10X, 'LOADING  CRITERION 7/10X ,' LOAD  =0  TRACTION  ONLY ’ /10X , ’ 

1L0AD  =1  ACCELERATION  ONLY '/10X, ’LOAD  =2  TRACTION  AND  ACCEL ' /10X , ' I 
2M  THIS  CASE  LOAD  =', I5//10X , 'LOADING  FUNCTION  IS '/10X , 'AZERO+BZERO 
3*C0S(WT)+CZER0*SIN(WT )  WHERE , '/10X , ’OMEGA= ' ,E15.4 , ' AZERO= ' ,E15 .4, ' 
4BZER0= ' ,E15.4,'CZER0=' ,E15.4//) 


FORMAT  (515) 

FORMAT  (10X,1 ACCELERATION  INPUT ' /10X , 1 LACFL  =0  LOADING  FUNCTION ‘ /I 
lOX.'LACEL  =1  ACCELEGRAM  RECORD  * /10X , 1 1 N  THIS  CASE  LACFL  = 1  ,15) 
FORMAT  (I0X, 'ACCELERATION  DEGREES  OF  FREEDOM') 

FORMAT  (10X.4I5) 

FORMAT  (215) 

FORMAT  ( // 1 0 X , ' TAPE  NUMBER  FOR  ACCELEGRAM  RECORDS- ‘ ,I5/10X, 'NUMBER 
1  OF  ACCELEGRAM  RECORDS  =',  1 5///10X  ,'  SUPPORT  ACCELERATION') 

FORMAT  (10X.E10.3) 

FORMAT  (4F10. 0,315) 

FORMAT  ( 3X , ' STEP  N0',2X,'END  TIME ', 5X ,' THETA ' ,6X ,' GAMA ' ,6X ,' BETA '  , 
IlX.'NO  OF  INC',  5X,  'PRINT', 4X , 'STRESS'/) 

FORMAT  (110, 4D 10.3,3110) 


FORMAT 

FORMAT 


'NUMBER 


'BETA'  , 


r-H  cnj  ro 
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SUBROUTINE  STEPS  (THETA, GAMA, BETA, DT ,R,F ,AK ,SK , ID ,C0EF ,MEQ , 

C  . . . . - 

1 IBAND ,NUMEL .NUMNP ,NDF,NST ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  AK(MEQ, IBAND),  R(NDF .NUMNP ) ,  ID(NDF,NUMNP ) ,  F(NDF,1 ) ,  SK 
1 (NST.NST) ,  SM(32 ,32) ,  SD(32,32),  LP(32),  C0EF(10) 

REWIND  16 
READ  (16)  AK,R 

C 

C  ASSEMBLE  THE  DYNAMIC  STIFFNESS  MATRIX 
C1=BETA*THETA*THETA*DT *DT 
C2=BETA*THETA*DT /GAMA 
REWIND  15 
DO  3  M=1 ,NUMEL 

READ  (15)  NSL ,  ( (SK{ I  ,J  )  ,J=1 ,NSL ) , 1=1 ,NSL ) , ( (SM( I ,J ) ,J=1 ,NSL ) , 1=1 ,N 
1SL ) , ( (SD ( I ,J ) ,J  =  1 ,NSL ) , I =1 ,NSL ) , (LP ( I ) , 1  =  1 ,NSL ) 

DO  2  1  =  1, NSL 
1 1 =LP  ( I ) 

DO  1  J=1,NSL 
JJ=LP ( J ) 

KK=JJ-I I +1 
IF(KK.LE.O)  GO  TO  1 

AK ( 1 1 ,KK )=AK (II ,KK )+SM(I ,J )/Cl+SD( I ,J )/C2 
CONTINUE 
CONTINUE 
CONTINUE 
C 

C  INTRODUCE  KINEMATIC  CONSTRAINTS 

DO  5  M=l, NUMNP 
DO  4  J  =  1  ,NDF 

IF(ID(J,M).EQ.O)  GO  TO  4 
IDF= (M-l )*NDF+J 

CALL  GEOMBC  (F(J ,M) , IDF, IBAND ,MEQ,R ,AK ) 

4  CONTINUE 

5  CONTINUE 
C 

C  OBTAIN  TIME  DOMAIN  INTEGRATION  FACTORS 

T1=THETA 
T2=T1*T1 
T3=T2*T1 
COEF ( 1  )  =  1 .D0/T3 
C0EF(2  )  =  1 . DO-COEF ( 1 ) 

C0EF(3 )=( 1 .DO-1 .D0/T2)*DT 

C0EF(4)=(1.D0-1.D0/T1)*UT*DT/2.D0 

C0EF(5)=GAMA/(BETA*DT*T3) 

C0EF(6 )=1 ,D0-GAMA/(BETA*T2 ) 

COEF (7 )= ( 1 ,D0-GAMA/(BETA*T1*2 .DO) )*DT 
COEF (8)=C0EF ( 5 )/(DT*GAMA) 

C0EF(9  )=COEF (5 )*T1/GAMA 

COEF ( 10 )  =  1 . DO— 1 .D0/(BETA*T1*2.D0) 

RETURN 

END 


SUBROUTINE  INCRS  ( THETA , GAMA .BETA ,DT ,R ,UO , VO ,AO,AN , BN , SK , 


1  ID, IX ,NDF, NUMNP .NUMEL ,NST .MEN ,NEN1 ,F, ABASE , OMEGA, AZERO  , 

2BZF.RO  ,C ZERO  .TIME  .LOAD  .LACE L  ,  I ACEL  .TRECD ,  ACEL  .NDTACL  ) 

IMPLICIT  REAL*8(A-H  ,0-Z) 

DIMENSION  R(NDF.NUMNP),  UO(NDF.l),  VO(NDF.l),  AO(NDF.l),  AN(NDF.l) 
1,  BM(NOF.l),  ID(NDF.NUMNP),  SK(NST.NST),  F (NDF.l ) ,  ABASE ( NDF , 1 ) ,  I 
2ACEL (4 ) ,  TRECD (20 ) ,  ACEL (20) ,  IX(NENl.l),  SM(32,32  ),  SD ( 32  ,32 ) ,  LP 
3(32) 

REWIND  17 
READ  (17)  R 

CALL  ZERO  (ABASE ,NDF*NUMNP) 

LOADING  FUNCTION  AT  TO+THETA*DT 

T=TIME+(THETA-1.D0)*DT 

ARGUM=OMEGA*T 

FACT =AZERO+BZERO*DCOS ( ARGUM )+CZERO*DS I N ( ARGUM ) 

IF (LOAD.EQ.O)  GO  TO  2 

ACCELERATION  LOADING 

IF(LACEL.NE.O)  CALL  RECORD  (T, TRECD, ACEL, NDTACL, ACLRCD) 

DO  1  1=1 .NUMNP 
DO  1  J=1,NDF 

IF(ID(J,I).NE.O)  GO  TO  1 
IF ( I ACE L (J).EQ.O)  GO  TO  1 
HARMONIC  ACCELERACTION 
IF(LACEL.EQ.O)  ABASE ( J , I )=FACT 
ACCELEGRAM  RECORDS  ACCELERACTION 
IF(LACEL.NE.O)  ABASE ( J  , I  )=ACLRCD 
CONTINUE 

I F ( LOAD . EQ . 1 )  GO  TO  4 
CONTINUE 

_ ADD  EXTERNALLY  APPLIED  CONC.  NODAL  LOADS  TO  R 

DO  3  N=l, NUMNP 
DO  3  J  =  1  ,NDF 

IF ( ID ( J ,N ) .NE .0 )  GO  TO  3 
R(J,N)=R(J,M)+FACT*F(J,N) 

CONTINUE 

CONTINUE 

X=THETA*DT 

XX=X*X 

Y=.5D0-BETA 

YY  =  1  .DO-BE  TA/GAMA 

ZZ= .5D0 -BETA /GAMA 

DO  5  1=1, NUMNP 

DO  5  J  =  1  ,NDF 

AN(J,I  )=(IJO(J,I  )+X*VO(J,I)+Y*XX*AO(J,I))/(BETA*XX) 

BN ( J , I )= (U0( J , I )+YY*X*VO (J , I )+ZZ*XX*AO ( J , I ) )*GAMA/( BETA*X ) 

REWIND  15 
DO  9  M=l, NUMEL 


READ  ( 15 )  NSL , ( ( SK ( I ,J ) ,J=1 ,NSL ) , 1=1 ,NSL ) , ( (SM( I ,J ) ,J=1 ,NSL ) , 1=1 ,N 
1SL ) , ( (SD ( I ,J ) , J  =  1 ,NSL ) , 1  =  1 ,NSL ) , (LP ( I ) , 1  =  1 ,NSL  ) 

NEL^NSL/NDF 
DO  8  I -1  ,NEL 
J = I X  ( I ,  M ) 

DO  7  I DF  =  1 ,NDF 

II=(I-1)*MDF+IDF 

IF ( ID( IDF ,J ) .NE .0 )  GO  TO  7 

DO  6  K=1,NEL 

N=IX(K,M) 

DO  6  L=1,NDF 
KK= ( K  —  1 ) *NDF+L 

R ( IDF , J )=R ( IDF , J )  +  SM ( 1 1 ,KK ) *( AN (L  ,N ) -ABASE (L ,N ) )  +  SD ( 1 1  ,KK  )*BN  (L  ,N  ) 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 


SUBROUTINE  INACL  ( AM ,U0 ,VO ,AO ,R ,F ,SK , IX , ID .MEN ,NEN1  ,NST  ,MEQ 

1 1  BAND , NUMNP , NUMEL , NDF , BZERO ) 

IMPLICIT  REAL*8(A-H,0-Z) 

LOGICAL  BACK ,FROW 

DIMENSION  AM(MEQ, IBAND ) ,  UO(NDF.l),  VO(NDF,l),  AO(NDF.l),  R(NQF, 1 ) 
1,  F(NDFJ),  SK(NST,NST),  SM(32,32),  SO(32  ,32  ),  LP(32),  IX(NENl.l), 
2  ID(NDF,1) 

THIS  SUBROUTINE  ASSEMPLE  THE  MASS  MATRIX  AND  SOLVE  THE 
EOIJATION  OF  MOTION  FOR  THE  INITIAL  ACCELERATION 
CALL  ZERO  (AM,IBAND*MEQ) 

CALL  ZERO  (AO ,NDF*NUMNP ) 

REWIND  15 
DO  6  M=1 , NUMEL 

READ  (15)  NSL ,  ( ( SK ( I  ,J )  ,J  =  1 ,NSL  ) , I  =  1 ,NSL  ) ,  ( ( SM( 1 ,0 ) , J  =  1 ,NSL ) , I  =  1 ,N 
1SL  ) ,  ( (SD ( I ,J ) ,J  =  1 ,NSL ) , 1  =  1 ,NSL ) , (LP ( I ) , I =1  ,NSL ) 

DO  2  1=1, NSL 
1 1 =LP ( I ) 

DO  1  J  =  1  ,NSL 
JJ=LP ( J ) 

KK  =J  J - 1  1+1 

IF(KK.LE.O)  GO  TO  1 

AM ( 1 1 ,KK )=AM( 1 1 ,KK )+SM(I ,J ) 

CONTINUE 

CONTINUE 

OBTAIN  THE  RIGHT  HAND  SIDE 
NEL=NSL/NDF 
DO  5  1  =  1  ,NEL 
J=IX( I ,M) 

DO  4  IDF=1 ,NDF 
1 1  =  ( 1-1 )*NDF+IDF 
DO  3  K  =  1  ,NEL 
N  =  I X ( K ,M) 

DO  3  L=1 ,NDF 
KK=(K-1)*NDF+L 

AO ( IDF ,  J  )  =  AO  ( IDF , J  )+SK  (II  ,KK  )*U0  (L  ,N  )+SD  ( 1 1  ,KK  )*V0 (L  ,N ) 

CONTINUE 

CONTINUE 

CONTINUE 


INTRODUCE  TRACTION  AT  TIME=0. 

DO  8  M=l, NUMNP 
DO  8  J=1,NDF 

IF ( I D ( J ,M).NE.O)  GO  TO  7 

AO(J ,M)=R(J ,M)-AO (J ,M)+3ZER0*F(J ,M) 

GO  TO  8 

AO  ( J  ,  M )  =  R  ( J  ,  M )  -  AO  ( J  ,  M ) 

CONTINUE 


INTRODUCE  CONSTRAINT  ACCELERATIONS 


ASUPR=O.ODO 
00  10  M=1 ,NUMNP 
00  9  J=1  ,NDF 

IF( ID(J.M).EQ.O)  GO  TO  9 
I  DF= ( M— 1 )*NDF+J 

CALL  GEOMBC  ( ASURP , IDF , I  BAND ,MEQ , AO , AM ) 
CONTINUE 


10  CONTINUE 


TRIANGULIZE  TH  MASS  MATRIX 
BACK=. FALSE. 

FROW=.TRUE . 

CALL  BANSOL  (AM ,A0 ,A0 ,MEQ, IBAMD .BACK ,FROW) 


SOLVE  FOR  INITIAL  ACCEL 
BACK=.TRUE. 

CALL  BANSOL  (AM,AO,AO,MEQ, IBAND.BACK ,FROW) 
RETURN 


SUBROUTINE  RECORD  (TIME .TRECD, ACEL .NDTACL .ACLRCD ) 


IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  TRECD ( 20 ) ,  ACEL ( 20 ) 

THIS  SUBROUTINE  INTERPOLATE  ACELOGRAM  DATA 

ACLRCD=O.ODO 

RATIO=O.ODO 

DO  1  I=2,NDTACL 

IF (TIME.GE. TRECD ( I ) )  GO  TO  2 

CONTINUE 

RETURN 

DTRCD=TRECD( I ) -TRECD ( 1-1 ) 

IF(DTRCD.EQ.O.ODO)  GO  TO  3 

RAT  1 0= ( TI ME— TRECD ( I  —  1 ) )/DTRCD 

ACLRCD=ACEL ( 1-1 )*( 1 .DO-RATIO )+ACEL ( I )*RATIO 

RETURN 


SUBROUTINE  ELMT02  ( IX ,XL ,D,FL ,S,UL ,TL ,NDM,NDF,NST,NEL , ISW 
1IL.LN  ,NG ) 


r fl| I  FD  RY 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /CDATA/  0 .HEAD (20 ) .NUMNP .NUMEL .NUMMAT ,NEN ,NEQ .MEQ.NTAPE ,NPN 
COMMON  /ANGLE/  ANG(4) ,DC(3) , ST ( 2 ) , IEPS .LINE.EWR 
COMMON  /GQDRT /  GPT(4 ,4) ,WT(4 ,4 ) ,SHP(4 ,9 ) ,LX(4 ) ,LE (4 ) 

DIMENSION  D( 1 ) ,  XL(NDM.l),  FL(NDF.l).  S(NST.NST),  UL(NDF.l),  TL ( 1 ) 
1,  I X ( 1 ) ,  SM (32,32) ,  SD(32,32) ,  LP(32) 


NGP=NG 

GO  TO  (1,1, 4, 5, 6),  ISW 

COMPUTE  ELEMNT  STIFFNESS,  MASS  AND  DAMPING  MATRICES 

CONTINUE 

NSL=NEL*NDF 

ELEMENT  LENGTH 

H=XL (NDM,2 )-XL ( NDM ,1 ) 

SOLID  STIFFNESS 

F=(D(1 )+D(19)*D(19)*D(18) )/H 

S(1 ,1  )=F 

S(3 ,3 )=F 

S(1 ,3)=-F 

SOLID  MASS 

F=D(23)*H/3.D0 

SM (1,1 )=F 

SM(3 ,3)=F 

SM(1 ,3)=F/2.D0 

COUPLING  STIFFNESS 

F=D(19 )*D(18 )/H 

S( 1 ,2  )=F 

S(  1 ,4 )=-F 

S(2,3)=-F 

S  ( 3 , 4  )=F 

COUPLING  MASS 

F=D ( 26 ) *H / ( D ( 2  2 ) *3 . DO ) 

SM ( 1 ,2 )=F 
SM ( 1 ,4)=F/2.D0 
SM(2,3)=F/2.D0 
SM { 3 , 4 ) =F 
FLUID  STIFFNESS 
F=D(18)/H 
S(2 ,2 }=F 
S(4,4)=F 
S  ( 2 ,4 ) =-F 
FLUID  MASS 

F=D(26  )*H/( D(22 )*D{22 )*3.D0) 

SM(2 ,2 )=F 
SM(4,4)=F 
SM(2 ,4)=F/2 
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00  2  I =1 ,NSL 
00  2  J  =  1  ,NSL 
S(J,I)=S(I,J) 

SM(J , I )=SM( I ,J ) 

2  SD( I ,J )=0.0D0 

C  FLU  10  0 AMP  INC, 

F=H/(D(24)*3.D0) 

SD(2 ,2  )=F 
SO (4 ,4 )=F 
SD(2,4)=F/2.D0 
SO (4 ,2 )=F/2.D0 
C  SOLID  DAMPING 

F 1 =0 (21) 

F2=D{22 )*D(22 )*D(21 ) 

F3=D { 20 ) 

F4=D(20)*D( 19  )*D ( 19  ) 

SD (1,1 )=F1*SM( 1 , 1 )-F2*SM(2 ,2 )+F3*S( 1 , 1 )-F4*S(2 ,2 ) 
SD(3,3)=F1*SM(3,3)-F2*SM(4,4)+F3*S(3,3)-F4*S(4,4) 

SD ( 1 ,3 )=F1*SM(1 ,3 )-F2*SM{2 ,4 )+F3*S( 1 ,3 )-F4*S(2 ,4 ) 

SD (3,1)  =  SD (1,3) 

DO  3  1=1, MEL 
11=2*1 

LP ( 1 1 ) = 2 * I X  ( I ) 

3  LP ( 1 1-1 )=LP ( 1 1  )-l 

WRITE  (15)  NSL,((S(I,J),J  =  1 ,NSL ) , 1  =  1 ,NSL ) , ( (SM( I ,J ) ,J  =  1 ,NSL  ) , 1=1  ,N 
1SL ) , ( (S0( I ,J  )  ,J  =  1 ,NSL  ) , 1  =  1  ,NSL ) ,  (LP ( I ) ,  1  =  1 ,NSL ) 

RETURN 

4  RETURN 

5  CONTINUE 

C  GET  FLUID  AND  SOLID  BODY  FORCES 

H=XL(NDM,2 ) — XL ( NDM , 1 ) 

RU0=D(23) 

RH0F=D(26) 

FL ( 1 ,1 )=H*RHO*DC (2 )/2 .DO 
FL (2 ,1 )=H*RHOF*UC (2 )/2 .DO 
FL ( 1 ,2  )=FL (1,1) 

FL ( 2  ,2)=FL (2 ,1 ) 

RETURN 

6  CONTINUE 

C  CALCULATE  EFFECTIVE  .TOTAL  STRESSES  AND  PORE  PRESSURES  AT 

C  ELEMENT  CENTRID 

H=XL (NDM.2 ) -XL (NDM, 1 ) 

XC=XL(NDM,1)+H/2.D0 
C  STRAIN  IN  X-DIR 

EPS=(UL(1,2)-UL(1,1))/H 
C  FLUID  VOL  STRAIN 

ZETA=(UL(2,2)-UL(2,1 ) )/H 
C  PORE  PRESSURE 

P=D(18)*(D(19)*EPS+ZETA) 

C  EFFECTIVE  STRESS  IN  X-DIR 

SIG=EPS*D( 1 )+ (D( 19 )-l .DO ) *P 
C  TOTAL  STRESS 


S I GT=S I G+D ( 19 )*P 
IF(LINE.NE.O)  GO  TO  7 
WRITE  (6,10)  HEAD 
PRINT  8 
LINE  =L I NE +1 
IF(LINE.EQ.50)  LINE=0 
PRINT  9,  LN ,XC ,SIG ,P ,SIGT 
RETURN 

FORMAT  (6X , 'EL ' ,8X , 'X-COORO ' ,6X , 'EFF-STRES ’ , 7X , ’POR-PRES ' ,6X , ‘TT 
1L-STRES 1 /) 

FORMAT  (5X.I5.4E15.4) 

FORMAT  ( 1H1 ,5X .20A4/5X , 1 _ ELEMENT  STRESS _ ’ , ‘ CCW  ROTATION  IS 

1  POSITIVE - * /5X , ‘ - 1-DIRECTION  IS  RADIAL  AND  2-DIRECTION  IS  AX 

2IAL  FOR  AN  1 ,1X , ‘AX  I  SYMMETRIC  PROBLEM _ '//) 


SUBROUTINE  ELMT03  ( IX ,XL ,D ,FL ,S ,UL ,TL ,NDM,NDF,NST ,NEL , ISW 
1 IL ,LN ,NGP  ) 


CALLED  BY  ELMLIB 

CALLS  SHAPE2 ,PVALUE 

IMPLICIT  REAL*8(A-H ,0-Z ) 

COMMON  /CDATA/  0 .HEAD (20 ) .NUMNP .NUMEL ,NUMMAT,NEN ,NEQ,MEQ,NTAPE ,NPN 
COMMON  /ANGLE/  ANG ( 4 ) ,DC { 3 ) ,ST (2 ) , IEPS .LINE , EWR 
COMMON  /GQDRT/  GPT (4 ,4 ) ,WT(4,4) ,SHP (4,9 ) ,LX (4 ) ,LE (4 ) 

COMMON  /STRES/  EPS (6 )  ,SIG(6 ) 

DIMENSION  D ( 1 ) ,  XL(NDM,1),  FL(NDF,1),  S(NST,NST),  UL(NDF,1),  TL ( 1 ) 
1,  IX(1),  SM (32,32),  SD(32,32),  LP (32 ) 

GO  TO  (1,1,8,9,11),  ISW 
L  CONTINUE 

....COMPUTE  ELEMENT  STIFFNESS,  MSASS  AND  DAMPING  MATRICES  BY (NGP*N 
NSL=NEL*NDF 
DO  2  1=1 ,NSL 
DO  2  J=1 ,NSL 
SD( I ,J )=0.0D0 
l  SM( I ,J )=0.0D0 

DO  5  M=1  ,NGP 
DO  5  N=1 ,NGP 

CALL  SHAPE2  (XL , IX ,NDM,NEL ,XDJ ,M,N ,NGP ) 

COLUMN 
Jl  =  l 

DO  4  J  =  1  ,NEL 
W1J=SHP (1 ,J )*XDJ 
W2J=SHP (2 ,J )*XDJ 
W3J=SHP(3,J )*XDJ 
K 1  =  J 1 
ROW 

DO  3  K=J  ,NEL 

S(J1  ,K1  )=S(J1 ,K1 )+WlJ*SHP(l ,K) 

S(J1,K1+1)=S(J1,K1+1 )+WlJ*SHP (2,K ) 

S(J1+1,K1)=S(J1+1,K1 )+W2J*SHP(l,K) 

S( Jl+1 ,K1+1 )=S(J1+1 ,K1+1 )+W2J*SHP(2,K) 

SM(J1 ,K1 )=SM(J1 ,K1 )+W3J*SHP(3,K ) 

3  K1=K1+N0F 

4  J 1 = J 1 +NDF 

5  CONTINUE 

D1=D(1 )+D(19)*D(19)*D(18) 

D2=D(2)+D(19)*D(19)*D(18) 

D3=D(3)+D(19)*D(19)*D(18) 

F1=D (21 ) 

F2=D (21 )*D{22 )*D(22 ) 

F3=D (20 ) 

F4=D ( 20 ) *D ( 1 9 ) *0 ( 1 9 ) 


DO  6  J=1,NSL,NDF 
DO  6  K=J  ,NSL  ,NDF 
W11=S(J ,K ) 

W12=S( J ,K+1 ) 

W21=S{ J+l ,K ) 

W22=S(J+1 ,K+1 ) 

A11=W11*D(18) 

A12=W12*D(18) 

A21=W21*D(18) 

A22=W22*D(18) 

B11=A11*D(19) 

B 1 2 = A 1 2  *D (19) 

B21=A21*D(19) 

B22=A22*D (19 ) 

DFF=SM{ J ,K ) 

WSS=DFF*D(23 ) 

WFF=DFF*D(26)/(D(22)*D(22) ) 

WSF =DFF*D ( 26 ) /D ( 22 ) 

SOLID  STIFFNESS 

S(J,K)=W11*D1+W22*D(4)+W12*D(5)+W21*D(5) 

S(J ,K+1 )=W12*D2+W21*D(4)+W11*D(5)+W22*D(6) 
S(J+1,K)=W21*D2+W12*D(4)+W22*D(6)+W11*D(5) 
S(J+1,K+1)=W22*D3+W11*D(4)+W21*D(6)+W12*D(6) 
S(K ,J )  =  S( J ,K ) 

S(K ,J+1 )  =  S( J  +  l  ,K ) 

S(K+1 ,J )  =  S(J ,K+1 ) 

S(K+1,J+1)=S(J+1 ,K+1 ) 

COUPLING  STIFFNESS 
S(J ,K+2 )  =  3 1 1 
S( J ,K+3 )=B12 
S( J+l ,K+2 )=B21 
S( J+l ,K+3 )=B22 
S ( J  +2 ,K  )=R11 
S(J+2,K+1  )=B12 
S ( J +3 ,K  )=B21 
S( J+3 ,K  +  1  )=B22 
S(K+2,J  )=B11 
S ( K +3 ,J  )=B12 
S ( K +2 ,J+1  )=B21 
S ( K+3 ,J+1  )=B22 
S(K ,J+2  )=B11 
S(K+1 ,J+2  )=B12 
S(K,J+3)=B21 
S ( K  + 1 ,J+3  )=B22 


FLUID  STIFFNESS 
S(J+2,K+2 ) = A 1 1 
S( J+2 ,K+3 )=A12 
S(J+3 ,K+2 )=A21 
S(J+3 ,K+3 )=A?2 


S(K+2  ,J+2 )=A1 1 
S(K+3  ,J+2 )=A12 
S(K+2  ,J+3  )=A21 
S(K+3 ,J+3 )=A22 

SOLID  MASS 
SM( J ,K )=WSS 
SM( J+l ,K+1 )=WSS 
SM(K ,J )=WSS 
SM(K+1,J+1)=WSS 

COUPLING  MASS 
SM (J+l ,K+3 )=WSF 
SM( J+2 ,K )=WSF 
SM( J+3 ,K+1 )=WSF 
SM(K+2 ,J )=WSF 
SM(K+3 ,J+1 )=WSF 
SM(K,J+2)=WSF 
SM(K+1 , J+3 )=WSF 
SM ( J ,K+2 )=WSF 

FLUID  MASS 
SM(J+2,K+2)=WFF 
SM( J+3 ,K+3 )=WFF 
SM(K+2 , J+2 )=WFF 
SM(K+3,J+3)=WFF 

FLUID  DAMPING 
SD(J+2,K+2)=DFF/D(24) 

SO (J+3 ,K+3)=DFF/D(25) 

SD(K+2 ,J+2 )=SD(J+2 ,K+2 ) 

SD(K+3,J+3)=SD(K+3,J+3) 

SOLID  DAMPING 

SD ( J  , K  )=F1*WSS-F2*WFF+F3*S( J ,K )-F4*S( J+2 ,K+2 ) 

SO(J ,K+1 )=F1*WSS-F2*WFF+F3*S(J ,K+1 )-F4*S( J+2 ,K+3 ) 

SD (J+l ,K )=F1*WSS-F2*WFF+F3*S(J+1 ,K )-F4*S(J+3 ,K+2 ) 

SD ( J+l ,K+1 )=F1*WSS-F2*WFF+F3*S(J+1 ,K+1 )-F4*S( J+3 ,K+3 ) 

SD(K,J)  =  SD(J,K) 

SD (K+l ,J )=SD(J ,K+1 ) 

SD(K ,J+1 )=SD (J+l ,K ) 

SD(K+1 , J+l )=SD (J+l ,K+1 ) 

DO  7  I =1 ,NEL 
1 1 =NDF*I 

LP (II )=NDF*IX( I ) 

LP ( 1 1-1 )=LP ( 1 1  )-l 
LP ( 1 1 -2 )=LP (II ) -2 
LP ( 1 1-3 )=LP ( II ) -3 

WRITE  (15)  NSL , ( (S( I ,J ) ,J  =  1 ,NSL ) , I =1 ,NSL ) , ( (SM( I ,J ) , J  =  1 ,NSL  ) ,  1=1  ,N 
1SL  ) , ( (SD( I ,J ) ,J  =  1 ,NSL ) , 1  =  1 ,NSL ) , (LP ( I ) , I =1 ,NSL ) 

RETURN 

CONTINUE 


....GET  FLUID  AND  SOLID  BODY  FORCES 
DO  10  M=1,NGP 
DO  10  N=1  ,NGP 

CALL  SHAPE2  (XL , IX ,NDM,HEL ,XDJ ,M,N ,NGP ) 

D2=XDJ*D(23 ) 

D4=XDJ *D(26 )/D(22 ) 

DO  10  1=1, NEL 
01=D2*SHP(3,I) 

D3=D4*SHP(3,I ) 

DO  10  J=1 ,2 

FL(J,I)=FL(J,I)+D1*DC(J) 

10  FL (J+2 , I )=FL ( J+2 , I )+D3*DC (J ) 

RETURN 

11  CONTINUE 

_ COMPUTE  STRESSES  AT  NGP  GAUSS  PTS. 

Z=O.ODO 
DO  14  M=1,NGP 
DO  14  N=1 ,NGP 

CALL  SHAPE2  (XL  , IX .NDM.NEL ,XSJ ,M,N ,NGP ) 

DO  12  1=1,3 

12  EPS ( I )=0.0D0 
ZETA=0.0D0 
X=0.0D0 
Y=0.0D0 

DO  13  J=1,NEL 
X=X+SHP(3,J)*XL(1,J) 

Y=Y+SHP (3 ,J )*XL(2 ,J ) 

EPS ( 1 )  =  SOLID  STRAIN  IN  X-DIRECTION  ,  EPSILON-XX 

EPS ( 2  )  =  SOLID  STRAIN  IN  XY -DIRECTION  ,  GAMMA-XY 

EPS ( 3 )  =  SOLID  STRAIN  IN  Y-DIRECTION  ,  EPSILON-YY 

EPSV  =  SOLID  VOLUMETRIC  STRAIN 
ZETA=  FLUID  VOLUMETRIC  STRAIN 
EPS(1)=EPS(1)+SHP(1,J)*UL(1,J) 

EPS(3 )=EPS(3  )+SHP (2 ,J  )*UL(2,J ) 

EPS (2 )=EPS(2 )+SHP (2 ,J )*UL { 1 , J )+SHP ( 1 ,J )*UL (2 ,  J ) 

13  ZETA=ZETA+SHP (1 ,J )*UL (3 ,J )+SHP (2 ,J )*UL (4 ,J ) 
EPSV=EPS(1 )+EPS(3) 

P  =PORE  PRESSURE 

P=D( 18 )*(D(19 )*EPSV+ZETA ) 

SIG(l)  =  STRESS  IN  X-DIRECTION  ,  SIGMA-XX 

S I G ( 2  )  =  STRESS  IN  XY-DIRECTION  ,  TAU-XY 

S I G ( 3 )  =  STRESS  IN  Y-DIRECTION  ,  SIGMA-YY 

....CALCULATE  EFFECTIVE  STRESSES 

SIG ( 1 )=D ( 1 )*EPS( 1 )+D(2 )*EPS(3)+(D(19 )-l .DO)*P 

SIG(2 )=D{4 )*EPS(2 } 

SIG(3)=D(2)*EPS(1)+D(3)*EPS(3)+(D(19)-1.D0)*P 
CALL  PVALUE  (SIG ,1 ,2 ,X ,Y ,Z  ,LN , IL,P ,D ) 


....CALCULATE  TOTAL  STRESSES 


FILMED 


